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ABSTRACT 


Numerical simulations using the thin-layer Navier-Stokes equations and chimera 
(overset) grid approach were carried out for flows around the integrated space shut- 
tle vehicle over a range of Mach numbers. Body-conforming grids were used for 
all the component grids. Testcases include a three-component overset grid — the 
external tank (ET), the solid rocket booster (SRB) and the orbiter (ORB), and a 
five-component overset grid — the ET, SRB, ORB, forward and aft attach hardware, 
configurations, The results were compared with the wind tunnel and flight data. 

In addition, a Poisson solution procedure (a special case of the vorticity-velocity 
formulation) using primitive variables was developed to solve three-dimensional, irro- 
tational, inviscid flows for single as well as overset grids. The solutions were validated 
by comparisons with other analytical or numerical solutions, and/or experimental re- 
sults for various geometries. The Poisson solution was also used as an initial guess 
for the thin-layer Navier-Stokes solution procedure to improve the efficiency of the 
numerical flow simulations. It was found that this approach resulted in roughly a 
30% CPU time savings as compared with the procedure solving the thin-layer Navier- 
Stokes equations from a uniform free stream flowfield. 
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1. INTRODUCTION 


1.1 Objective 


Over the years, many researchers have been searching for ways to simulate flows 
around complex geometries, such as full aircraft configurations (Jameson et al., 1986; 
Flores et al., 1987), and the integrated space shuttle configuration (Szema et al., 
1988; Buning et al., 1988, 1989). This goal demands the combined use of different 
technologies developed for computational fluid dynamics and other fields, like grid 
generation, flow solver, and flow visualization through the use of dedicated graphics 
workstations. The final results must be capable of capturing the significant flow fea- 
tures, and be verified with experimental data. This report describes the computation 
of the flow around the launch configuration of the integrated space shuttle vehicle 
as well as ways to improve the efficiency of numerical methods for flow simulations. 
The primary components of the configuration consist of the external tank (ET), the 
solid rocket booster (SRB), and the orbiter (ORB); and the secondary parts include 
the forward and aft attach hardware which link the external tank and the orbiter. 
The computed solutions were generally in good agreement with the wind tunnel tests 
and the flight data. With this achieved, the flow simulation around the space shuttle 
geometry can be further refined to include the parts purposely excluded in the cur- 
rent research, for example, the ORB vertical tail, the ET/SRB attach ring, and the 
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like. These were excluded because of limitations on modeling capability and compu- 
tational resources. With increasingly faster computers and continuous development 
of computational fluid dynamics, it will not be long before the quasi-static numerical 
simulation of the flow around complex geometries through different flow regimes (or 
Mach numbers) becomes a routine job for the major aircraft manufacturers and is 
integrated into the design loop. Hopefully, this technology will someday advance to a 
level that real-time simulation becomes a reality and the information rendered by the 
computer can help in real-time control on increasingly complicated flying machines. 

The primary objective of computer simulations of the space shuttle is to supple- 
ment the available experimental and flight data which suffer from inadequacies due to 
scaling effects, wind tunnel wall-interference effects, sting interference effects, instru- 
mentation limitations and the difficulty of safely obtaining valid flight data. Since the 
computer simulations are quite flexible, in that they allow for easy reconfiguration of 
the shuttle geometry, subject to the limitations of the flow solver and the difficulties 
of gndding, the simulations can easily be carried out for different geometries and 
different flow conditions. In this study, the simulation had been carried out from 
a grid consisting only of the ET, SRB and ORB to a more refined model with the 
addition of the forward and aft attach hardware, and Mach numbers ranging from 0.6 
to 2.0. The results from the simulations have already helped in the diagnostics of the 
damaged thermal protection system on the ORB surface on one of the shuttle mis- 
sions (Li, 1989). Li reported using PLOT3D (Buning and Steger, 1985), a graphics 
tool developed by Buning of NASA Ames Research Center, to draw particle traces 
based on the computed solutions from the damaged area to determine the possible 
debris path during the launching period. Other possible applications from the nu- 
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merical simulations are to predict the aerodynamic behavior during the emergency 
abort maneuvers or to determine the escape path for the astronauts. All of these 
applications are either impossible, impractical or extremely expensive to evaluate in 
the wind tunnel or in real flights; while for computer simulations, once the solution 
is obtained, the complete three-dimensional flowfield can be analyzed with the aid 
of modern graphics software and hardware to provide additional insight that might 
have otherwise been neglected, difficult, or impossible to capture in an experimental 
setup. The computer simulation can also be used as an aid to check the validity of 
the aerodynamic data base for the space shuttle; for example, the computed wing 
load from the experimental data may be different from what was observed in the real 
flight; thus, the numerical solution provides a third check for the data base. In the 
design or modification of the shuttle geometry, the relatively low cost computer sim- 
ulation (as compared with the cost of the wind tunnel test) allows room for designers 
to conduct various numerical experiments with a reasonable cost as is currently being 
done at Rockwell International for the ET/SRB attach ring. 

1.2 Approach 

This research was carried out as part of the space shuttle flow simulation project 
in the Applied Computational Fluids Branch (RFA) at NASA Ames Research Center. 
In one phase of the work, the thin-layer Navier-Stokes equations were solved using 
an implicit approximately factored finite-difference procedure (Steger et al., 1986; 
Ying et al., 1986) for the flows around all the components of the integrated space 
vehicle during its ascent mode for various nominal and abort flight conditions. The 
Baldwin-Lomax turbulence model was used to calculate the turbulent eddy viscosity. 
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Due to the complexity of the geometry, the chimera (Steger et al., 1983; Benek et al., 
1985, 1987; Dougherty et al., 1985) composite-grid approach was chosen to manage 
data communication between different grids. For each component (the external tank, 
the solid rocket booster, the orbiter and the forward and aft attach hardware) the 
body-conforming grid was generated independently and overlaid with each other to 
form the composite grid for the shuttle launch configuration. More details on the 
grid will be given in the chapter on the chimera approach. 

Since the geometry of the entire space vehicle is very complicated, it is extremely 
difficult, if not impossible, to model every possible detail of the entire vehicle. Thus, 
some simplification was made on the geometry, such as the elimination or idealization 
of the ET/ORB forward and aft attach hardware, the elimination of the ET/SRB 
attach ring, ET sway braces, the main fuel and oxidizer feed lines, etc. In addition, 
certain engineering approximations were made. For instance, stings have been used 
to represent plumes. This is a reasonable approximation for supersonic free stream 
flow due to the limited upstream influence; while at lower Mach numbers, the shuttle 
is in the lower, more dense atmosphere, and the plume expansions tend to be greatly 
reduced. 

In the first stage, a coarse grid (about 250,000 points for the whole shuttle 
configuration) was used for a flow at a free stream Mach number of 2.0. Since this 
was a supersonic flow, the required computer time was relatively small as compared 
to a subsonic flow due to the lack of upstream signal propagation, and was a proper 
testcase for testing the feasibility of the entire numerical procedure. 

In the second stage, some modifications were made in the geometry. For example, 
the ET/ORB attach hardware was added with some idealization of the geometry, 



5 


and the ET sting was removed. Various Mach numbers and different angles of elevon 
deflection for the ORB were computed and compared with the experimental and 
flight data. All the calculations in the first two stages were carried out on the NASA 
Numerical Aerodynamic Simulation (NAS) Cray 2 computer. 

Though it is possible to simulate flows around complex geometries such as the 
one presented in this research, the required computer time (6 to 20 hours on the Cray 
2 computer for the shuttle geometry) is still inhibitively high even for large aircraft 
manufacturers. This prompted the research on a faster procedure for the numerical 
simulation over complex geometries. At this point, a vorticity-velocity formulation 
was implemented to obtain a rough estimate of the flowfield 1 and the solution was 
then fed back to a more accurate solver for further calculations. The reason that this 
procedure may reduce the overall computing effort is that most flow solvers utilize 
a great amount of computer time in settling out the transient state; therefore, if 
a rough estimate of the flowfield is available, it is possible to reduce the computer 
time in obtaining the flow solution. The other possible use of the formulation is 
to evaluate a grid without putting in great effort just to find out that the grid is 
not adequate for the intended flow geometry. The solution of this formulation was 
carried out for several generic geometries (sphere, ellipsoid, and the external tank) 
to verify the validity of the formulation. Finally a solution for the entire shuttle 
vehicle was presented. All solutions from this formulation were compared with either 
the analytical solution or the solution obtained from the thin-layer Navier-Stokes 
equations as mentioned above. 


^Only inviscid flows were solved using the vorticity-velocity formulation. 
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1.3 Literature Review 

Navier-Stokes equations or its simplified approximations, such as the Euler equa- 
tions which neglect the viscous terms, have been used to solve for flows in many 
engineering applications. Depending on the choice of the dependent variables, the 
mathematical formulations of the Navier-Stokes equations (or its simplified approxi- 
mations) may be divided into the following categories (Guj and Stella, 1988), 

• vorticity/stream-function (Fromm and Harlow, 1963; Benjamin and Denny, 
1979), which in three-dimensional flows extends to vorticity-vector potential 
(Mallinson and de Vahl Davis, 1973; Richardson and Cornish, 1977) 

• vorticity / velocity 

• primitive variables. 

The focus here will be on the formulations suitable for solving three-dimensional flows. 
Thus, the review on solving the two-dimensional flows using the vorticity/stream- 
function method is left out; the reader is referred to Roache (1972) for a good review 
on the method using this formulation. 

1.3.1 Vorticity/velocity-vector potential 

The vorticity/stream-function formulation has been one of the popular approaches 
used for two-dimensional flow simulations. While the three-dimensional counterpart 
of the stream-function, referred to as the vector potential, or velocity-vector poten- 
tial, has been known for over a century (Helmholtz’s decomposition theorem, 1858, 
1867), it was not successfully implemented in computations until 1967 (Aziz and Hel- 
iums, 1967). The reasons for such a long delay were the difficulty in finding the proper 
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boundary conditions for the vector potential and the demands on computer resources 
since this formulation requires solutions for more variables than the primitive vari- 
able formulation. Aziz and Heliums (1967) in their implementation utilized boundary 
conditions similar to those formulated for a general hydrodynamic flowfield by Hi- 
rasaki and Heliums (1968) to study three-dimensional laminar natural convection in 
enclosures. The study of Aziz and Heliums showed that the vorticity/vector poten- 
tial approach can lead to faster and more stable convergence than for the comparable 
primitive variable formulation. Although the boundary conditions presented in Hi- 
rasaki and Heliums (1968) were fairly complete, their complexity in treating through- 
flows rendered the method useless in these situations. Later, Hirasaki and Heliums 
(1970) proposed a simplified boundary condition with the introduction of a scalar 
potential to account for the through-flow velocities. This formulation, termed “dual 
potential” (Chaderjian and Steger 1983), has been used to solve three-dimensional 
natural convection in enclosures (Mallinson and de Vahl Davis, 1973; Ozoe et al., 
1976, 1977, 1979, 1985) and three-dimensional flows in ducts (Aregbesola and Bur- 
ley, 1977; Wong and Reizes, 1984, 1986). 

Due to the added variable, i.e., the scalar potential, in Hirasaki and Heliums 
(1970), the vorticity/vector potential approach suffered from extra demands on com- 
puter time and storage. Wong and Reizes (1984) used the normal component of the 
specified inlet velocity vector in place of the scalar potential to reduce the requirement 
of computer storage. Calculations using this formulation for flows in a constant cross 
section duct were performed for a wide range of Reynolds numbers. Later, Wong and 
Reizes (1986) extended this approach to flows in multiply connected regions. 

Efforts have also been made to extend the vorticity/vector potential formula- 
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tion to deal with compressible flows. Hafez and Lovell (1981, 1983, 1988) devised an 
entropy and vorticity correction procedure for the potential/stream-function formula- 
tion. They showed details of the treatment of shocks and wakes, and the solution was 
compared with Euler solutions. Rao et al. (1987, 1989) combined the boundary-layer 
equations with the vorticity /vector potential formulation to do viscous-inviscid in- 
teraction. The vorticity was injected from the boundary-layer edge into the potential 
flow region to obtain the viscous effect. Gegg (1989) extended the vorticity /vector 
potential method to solve through-flow problems with heat transfer. 

1.3.2 Vorticity/velocity 

Only a handful of research projects have been carried out using the vortic- 
ity/velocity formulation for flow simulations. The earliest calculation using this for- 
mulation was reported in Fasel (1976). At each time step, Fasel solved two Poisson 
equations, derivable from the definition of vorticity, for the components of the ve- 
locity vector for a two-dimensional flow. This approach was applied to the study of 
the stability of boundary-layers. Dennis et al. (1979) extended the two-dimensional 
vorticity/velocity approach to three-dimensional steady flows by solving the Navier- 
Stokes equations in a cubical driven box. His approach bore some similarity to that 
of Aziz and Heliums (1967). The main difference was that three equations connecting 
the velocity and vorticity components were employed instead of the vector potential. 
These three equations, together with the three vorticity transport equations, form six 
simultaneous second-order partial differential equations to be solved. Gatski et al. 
(1982) applied compact finite-difference schemes to the vorticity- velocity form of the 
two-dimensional Navier-Stokes equations. The numerical solutions were obtained for 
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driven cavity flows. Fasel and Booz (1984) investigated the axisymmetric supercritical 
Taylor vortex flow for a wide gap. Farouk and Fusegi (1985) studied the natural and 
forced convection and heat transfer in a two-dimensional annulus. A coupled solution 
procedure was used for solving simultaneously the dependent variables using a block 
tridiagonal matrix inversion algorithm. The formulation was found to be fairly stable 
over a large range of Reynolds and Rayleigh numbers. Orlandi (1987) solved high- Re 
flows using a block ADI method which strongly coupled field equations and boundary 
conditions and satisfied the continuity equation without requiring an iterative proce- 
dure. Driven cavity and backward facing step flows were computed and verified with 
other numerical solutions and/or experimental results. Osswald et al. (1987) solved 
the three-dimensional unsteady Navier-Stokes equations using a direct inversion pro- 
cedure for a shear driven viscous flow within a cubical box. The three-dimensional 
vorticity transport equation in their procedure was solved using an approximate fac- 
torization method which required the inversion of only scalar tridiagonal matrices, 
rather than the usual block-tridiagonal systems. Guj and Stella (1988) computed 
two-dimensional incompressible flows in driven cavity and over a backward-facing 
step using a scalar ADI method. Speziale (1987) elaborated on the advantage of the 
vorticity /velocity formulation in a non-inertia coordinate system. He showed that 
the non-inertia effects, arising from both the rotation and translation of the frame 
of reference relative to an inertia frame, only enter into the equation through the 
implementation of initial and boundary conditions. This is in contrast to the primi- 
tive variable formulation, where non-inertia effects appear directly in the momentum 
equations in the form of Coriolis and Eulerian accelerations which may give rise to a 
variety of numerical problems (Williams, 1969). Considering the relatively short list 
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of references for this formulation, it is evident that the vorticity/ velocity approach 
needs further study. The major reason for the lack of research on this formulation 
may be attributed to the extra variables needed which in turn translate into extra 
demand on computer storage. For three-dimensional flows, memory requirements 
tend to govern the feasibility of implementing a particular numerical scheme. It is 
therefore understandable that numerical schemes which deal with a smaller number 
of variables have been strongly favored. 

Of Jill the studies mentioned above for the vorticity /velocity formulation, none 
were conducted for compressible flows. In the current study, the Crocco relationship 
was used in place of the vorticity transport equations; however, it is valid only for 
inviscid flows. Bernoulli’s equation for compressible, non-isentropic flows was used 
to account for density changes. Although only irrotational flow was computed in the 
current research, the proposed scheme has the potential of treating inviscid rotational 
flows. For the computation of viscous flows, the vorticity /velocity procedure imple- 
mented in the current research can be coupled with the boundary-layer (as in Rao 
et al., 1987, 1989) or Navier-Stokes equations to solve the viscous flow. The viscous 
effects in the inviscid region are accounted for through the injection of vorticity from 
the boundary-layer edge or the computational domain of the Navier-Stokes equations. 

1.3.3 Primitive variables 

Numerous methods using primitive variables have been developed to solve the 
Navier-Stokes equations. However, only methods dealing with improving the effi- 
ciency of the numerical flow simulation and techniques to treat complex geometries 
will be given attention here. Reviews of other related numerical methods for solving 
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the Navier-Stokes equations can be found in Shang (1985) and Holst (1987). An 
outlook for computational aerodynamics was given by Chapman (1979). 

With the debut of increasingly powerful computers, the quest for numerical flow 
simulations around real-world geometries, e.g., a complete aircraft, is within reach for 
those with access to supercomputers. Jameson and Baker (1987) computed an Euler 
solution around a Boeing 747-200 using an unstructured grid consisting of tetrahedral 
meshes. Obayashi (1987) carried out a Navier-Stokes solution for the ONERA M-5 
model (a wing- fuselage-tail geometry) using a single grid system. Flores et al. (1987) 
reported a zonal approach to obtain a transonic flow solution to the Navier-Stokes 
equations for a fighter-like configuration, while Buning et al. (1988) used an overset 
grid approach to solve the thin-layer Navier-Stokes equations for the integrated space 
vehicle launching configuration over a range of Mach numbers. Of those reported 
calculations for complex geometries, the grid systems used can roughly be divided 
into the following types: 

• single grid 

• multiple grids — can further be divided into 

1. zonal (patched) grid 

2. overset grid 

• adaptive grid. 

Often, a single grid generated for a complex geometry contains overly skewed 
meshes which in turn give rise to inaccurate solutions. Thus, a significant amount 
of effort is usually needed to modify the existing grid generater to yield acceptable 
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grids. In addition, the single grid approach does not have the flexibility of selective 
grid refinement and may require more points to resolve the flow than a comparable 
multiple-grid method. While the multiple grid approach is flexible, it is not without 
its own problems. For the patched grid, each subgrid is generated subject to boundary 
constraints placed by the neighboring grids, and multiple grids must be interfaced 
and managed (Steger and Benek, 1986). Rai (1986a, 1986b) developed a scheme 
such that the zonal boundaries were treated in a conservative manner so that the 
discontinuities could move freely across these boundaries. Many calculations were 
carried out using patched grids to resolve gradients, treat moving boundaries (Rai, 
1985), and complex geometries (Eberle and Misegades, 1986). 

For the overset grids, the grid does not need common boundaries between sub- 
grids, but rather, a common or overlap region is required to allow ways for matching 
the solutions across boundary interfaces. Usually, interpolation is used for the so- 
lution matching among subgrids; however, this will not ensure conservation of flux 
quantities, and inaccuracies can occur in shock capturing. Each subgrid in this ap- 
proach is generated independently, which, in turn, reduces a complex grid generation 
problem into a series of simple ones. Atta and Vadyak (1982) devised this approach to 
solve the potential equation for two- and three-dimensional flows. Steger et al. (1983) 
and Benek et al. (1983) developed a chimera scheme to solve the two-dimensional 
Euler equations. Subsequently, the scheme was extended to treat three-dimensional 
flows (Benek et al., 1985). Buning et al. (1988) used the scheme to solve for flows 
around the integrated space shuttle vehicle. A related application by Wedan and 
South (1983) employed a Cartesian mesh in which the body was embedded. 

The adaptive grid method allows the mesh to evolve with the solutions and does 
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not need, for the initial mesh, to anticipate accurately the large gradient regions. An 
advantage of this approach is that the grid points are efficiently used and human in- 
tervention is not needed to place the grid points in regions of large gradients. Gnoffo 
(1982) modeled the mesh as a network of springs with the spring constants deter- 
mined by the gradient of flow variables. Ghia et al. (1983) coupled the grid-evolution 
equation to the flow equation by requiring that the coefficients of the convective 
terms be minimized. Brackbill (1982) and Saltzman and Brackbill (1982) used a 
variational technique to produce grid-evolution equations. Berger and Oliger (1984) 
developed a dynamic refinement method which embeds finer and finer grids to resolve 
flow gradients. 

Although it is feasible to carry out flow simulations on the current generation 
of supercomputers for complex geometries like a complete aircraft, the work still 
demands a significant amount of computer resources. Therefore, it is still necessary 
to improve the rate of convergence for the numerical algorithms. Van Dalsem and 
Steger (1985) implemented a “fortified” Navier-Stokes approach, in which solutions 
to the subset equations, e.g., the boundary-layer equations, were used to add forcing 
terms to the Navier-Stokes algorithm in the proper flow regions. This approach 
was found to improve the efficiency as well as the accuracy of a given Navier-Stokes 
algorithm. Van Dalsem and Steger (1986) solved the boundary-layer equations on a 
fine grid near the wall to resolve the viscous gradients near the wall. The boundary- 
layer solution was then used as a forcing function and interpolated to the coarse grid 
solved by a Navier-Stokes algorithm. They reported a 20-fold increase for the rate of 
convergence in their testcases. 

Another approach is to use different sets of equations for different flow regions, 
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e.g., the boundary-layer equations for the boundary-layer region and the potential 
equation for the potential flow far from the body. Various combinations of equa- 
tions are possible. Whitfield et al. (1981) used the Euler equations to obtain the 
inviscid flow solution and boundary-layer equations in the viscous layer. Halim and 
Hafez (1984) developed a scheme in which the stream function was used to calculate 
the inviscid flow and the partially parabolized Navier-Stokes equations were used for 
the near wall shear layer. The third and the most commonly used approach is the 
boundary-layer equations for the viscous layer and the full-potential equations for 
the inviscid flow. The two solutions are matched by iterating for the displacement 
thickness. Most of these have concentrated on solving the two sets of equations simul- 
taneously to improve the rate of convergence (Lee and Pletcher, 1986). Although the 
approach has taken into account the physics at different flow regions and can possibly 
save a significant amount of computer time if each flow region is resolved properly, 
the complexity of treating multiple solution algorithms and domain interfaces usually 
makes the coding more difficult. 
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2. CHIMERA APPROACH AND GRID TOPOLOGY 


2.1 Introduction 


For complex geometries, generating the grid for the flow solver is itself a difficult 
task. Though it is possible to generate a single grid for a complex geometry, the 
resultant grid is most often overly skewed in one direction or another, or doesn’t 
have the needed clustering to resolve the flowfield in regions of rapid change. A 
natural way to overcome this difficulty is to divide the complex shape into several 
simple ones and generate the grid about these simple shapes, then either patch them 
together, the so-called patch grid, or overlay one on top of the others, as an overset 
grid. Combined with inter-grid communication in the flow solver, the patched grid 
or overset grid can be used to treat complex geometries. Since this research was done 
entirely with overset grids, the details of the gridding are described only for the overset 
grid. The overset grid approach used in this research was first devised by Steger et al. 
(1983) and given the name, the chimera approach, after the Greek legendary creature 
that was compounded of incompatible parts which signifies that the chimera approach 
can take incompatible grids (i.e., no common boundaries between different grids) and 
“glue” them together to be solved by the flow solver. 

Although drawbacks exist in the composite grid, patched or overset, like difficul- 
ties in accurately passing boundary data in-between sub-grids as well as finding the 
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interfacing information at the boundaries that separate different component grids, 
several advantages of applying the chimera approach outweigh the concerns of the 
drawbacks mentioned. First, the chimera approach does not require common bound- 
aries between component meshes. Due to this characteristic, the component grid can 
be generated separately, thus degenerating the complexity of grid generation into a 
combination of a series of simple ones. Furthermore, changes of some component 
grids do not usually involve changes of other component grids, thus allowing more 
flexibility than other approaches in constructing a grid for a complex geometry. This 
approach also saves the effort of gridding for a complex geometry since the chimera 
approach allows for arbitrarily adding or subtracting component grids. For instance, 
the gridding for the complex geometry can start with a simple one for initial testing 
and gradually be refined to a more accurate one by adding more component grids 
without regenerating from scratch the whole grid for the complete configuration. Be- 
sides, each component grid can be tested individually and added later when it is 
good enough or needed. The second feature of the chimera approach is that the flow 
simulation is done in sequence for component grids. This approach offers a savings 
in memory usage for solving flows around a complex geometry since it only requires 
memory enough to handle the largest component grid. The nature of solving each 
component grid in sequence also suggests the possibility of using different schemes, 
e.g., different set of equations, different time steps, etc., for different components. 
This opens many possibilities for enhancing the rate of convergence. For example, in 
the case of steady flow, it is possible to carry out more iterations for the component 
grid having the slower rate of convergence or using different time steps for different 
component grids. The chimera approach is also readily available for multitasking if 
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sizes of the component grids are roughly the same due to its “separate and conquer 
approach. Or when the complex geometry consists of one large grid and several 
much smaller grids, it is possible to carry out the computation for the large grid on 
one processor and the rest of the smaller grids in sequence on another processor or 

processors. 


2.2 Chimera Approach 

Several important concepts and implementation details underlying the chimera 
scheme are briefly described in this section. For more details the reader is referred 
to Steger et al. (1983), and Benek et al. (1983, 1985-1987). The chimera scheme 
involves the composite of the overlapping grids (generated individually), and the 
intergrid communications. As each mesh is generated individually, some grid points 
in one mesh will inevitably fall within the body boundary of another grid or grids, 
thus creating one or more “holes” in the mesh. The ET grid shown in Figure 2.8a 
depicts the holes created by the presence of the ORB and the SRB; the points in the 
ET grid surrounding the ORB and the SRB are called the hole boundary points. The 
values of the flow variables at these points are interpolated from solutions on either 
the ORB or the SRB grid, thus creating a link between the ET grid and grids of the 
ORB or the SRB. All the points within the hole boundary (including the boundary 
itself) will not enter into the solution process via a flag (will explain later) in the flow 
solver to differentiate them from the field points. In Figures 2.8b and 2.8c, not only 
the hole boundaries but also the outer boundaries of the ORB and the SRB grids are 
used to establish links to other grids. The composite of the overlapped grids and the 
interpolation data at the interface boundaries among component grids are created by 
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a code, named PEGASUS, developed at CALSPAN of AEDC (Arnold Engineering 
Development Center). 

2.2.1 Hole creation 

As explained above, the holes in a mesh are due to the presence of the solid 
bodies embedded in the mesh. To find the hole points as well as the hole boundary, the 
current implementation basically involves a two-step procedure. First, by introducing 
an imaginary rectangular box enclosing the embedded body, all the points that fall 
outside of the sphere with diameter equal to the diagonal of this rectangular box are 
considered field points. This method is fast, though somewhat crude, and cheap as 
compared with the method (will be explained below) used to find the hole points. 
Thus, it is used to filter out most of the points from the hole searching procedure. 
If the points tested fall within the sphere, they may be inside the embedded body. 
A more accurate method is needed to tell whether or not a point is a hole point. 
To clarify the basic idea underlying the hole searching procedure, a two-dimensional 
instead of a three-dimensional case is presented to avoid unnecessary confusion. In 
Figure 2.1, after the point, P, is tested and found to fall within the sphere mentioned 
above, the point nearest to P is found on the surface of the embedded body, say 
P c , and from this point, an outward normal, N, is constructed. If the dot product, 
N ■ Rp < 0, P lies within the hole and a flag variable, IBLANK, is set to zero; 
otherwise, P is outside the hole and IBLANK is set to 1. The IBLANK variable is 
used by the flow solver to determine whether a point should enter into the solution 
process or not as illustrated by the equation below. 


AAQ = IBLANK * RHS 


( 2 . 1 ) 
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where A is the coefficient matrix, AQ the change in flow solution, and RHS the source 
term. Thus, for the hole points, the above equation reduces to 

AQ = 0 ( 2 - 2 ) 

and the values of the variables at the hole points are not changed in the solution 
process. After the hole points are found, the hole boundary points can easily be 
located by searching the IBLANK values of neighboring points. If any neighboring 
points have a zero IBLANK value, they are defined as hole boundary points and 
also assigned zero as their IBLANK values since they are updated from the embed- 
ded mesh and should not enter into the solution process. Figure 2.2 illustrates the 
searching procedure for the hole boundary points. Note that the updating procedure 
for the hole boundary and the outer boundary is explicit and may somewhat affect 
the convergence and stability of the numerical scheme. 

2.2.2 Interpolation points 

Values of the flow variables at the hole boundary points are interpolated from 
the embedded mesh. Thus, it is necessary to find the interpolation points on the 
embedded mesh from which the hole boundary values are interpolated. The procedure 
involves locating the point in the embedded mesh closest to each hole boundary point, 
and once the closest point is found, its pointers (array indices) are added to a list of 
such points to be used by the flow solver to update the variables at the hole boundary 
points. To reduce the effort in finding the closest point, each search is started from 
the point found to be the nearest in the previous search. The interpolation points 
for the outer boundary are found with the same procedure. 
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Figure 2.1: 


Method for locating points within a hole 



Figure 2.2: Hole boundary point construction 
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2.2.3 Interpolation 

The transfer of information for the overset grid is through the hole and outer 
boundaries. The accuracy of the interpolation procedure influences the accuracy of 
the solution. Steger et al. (1983) reported a significant amount of mismatch in the 
vicinity of the shock/grid boundary intersection for an airfoil using single and two- 
• grid configurations. The interpolation scheme, which was based on a Taylor series 
expansion, was suspected to be the cause. Mastin and McConnaughey (1984) showed 
that bilinear interpolation in two dimensions is better than Taylor series expansions 
when higher order derivatives of the solution are not important. In the current 
research, the trilinear interpolation of the form: 

<f> = ag + + a2V + a 3^ + a 4 Zl + a 5^C + a QVC + a 7^vC (2.3) 

was used. In the above equation, 0 < £,r/,< < 1 are the coordinates of the point to 
be interpolated and ag to a 7 are computed based on the values at the points forming 
the interpolation stencil. 


2.3 Grid Generation 

Since the real geometry of the integrated vehicle is very complicated, it is im- 
possible to include all the details in the computer flow simulation with the current 
state of technology and limited computer resources. As evidenced in Figure 2.3, the 
ET fuel feed lines, the ET/SRB attach ring, the ORB vertical tails, and the space 
shuttle main engine (SSME) are all clearly visible and may influence the surrounding 
flow. However, to demonstrate the feasibility of the numerical model, the geometry 
of the integrated space shuttle vehicle was simplified and idealized to a certain degree 
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in the early stage of the research. Specifically, all three major components, ET, SRB, 
and ORB, were modeled with stings at the back extending to the outflow boundary 
and the ORB without the vertical tail as shown in Figure 2.4. The first calculation 
was for a free stream Mach number of 2.0. For this flow, the upstream influence was 
small, so that it was possible to capture meaningful flow phenomena despite the geo- 
metric simplifications. Later, the attach hardware, forward and aft, were added and 
the ET sting was removed to more accurately model the real geometry as illustrated 
in Figure 2.5. The stings behind the SRB and the ORB were still kept to mimic the 
effect of the plume. The elevons of the ORB were also deflected to the wind tunnel 
testing or real flight position. The surface definition of the integrated shuttle vehicle 
was provided by Ben-Shmuel of Rockwell International. 

Body-conforming grids were used for all the component grids and the grid lines 
were clustered near the body surface to resolve the high gradients in the boundary 
layer. Not only does the use of body-conforming grids make the boundary condi- 
tions simple to implement, but it also facilitates the clustering of the grid points in 
the boundary layer. In general, the grid is mapped onto a uniformly spaced com- 
putational domain, (£, 77 , £), with £ aligned with the major flow direction, 77 in the 
circumferential direction, and ( away from the body. The orbiter grid is generated 
using a three-dimensional hyperbolic grid generator developed by Steger and Rizk 
(1985). The hyperbolic grid generator basically solves three simultaneous partial dif- 
ferential equations — two orthogonality relations between £ and £ and between 77 and 

C» 

re - tv = 0 (2.4) 


F77 • tv = 0 


(2.5) 
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and a user specified volume constraint, 


d{x,y,z) 

d{i,Th() 


= AV" 


( 2 . 6 ) 


where f is the position vector, (x,t/,z) f . The grid is obtained by first defining a 
surface grid, then using a hyperbolic grid generator which marches in the outward 
normal direction from the given surface distribution through the constraint of the 
two orthogonal relations and the user specified spacings (volume). Due to the nature 
of the time like marching in the outward normal direction, the location of the outer 
boundary can not be specified. However, for external flows, the location of the outer 


boundary does not have to be fixed at a predetermined location as it does for internal 


flows. Details of the hyperbolic grid generation procedure can be found in Steger and 
Rizk (1985), while the specific details related to the ORB grid generation are given 
in Rizk and Ben-Shmuel (1985) and Rizk et al. (1985). Figure 2.6 shows the different 


views of the ORB grid. 

For the ET and SRB, the grids were generated using a two-dimensional hy- 
perbolic grid generator since the geometries were axisymmetric and could be spun 
around 360° to obtain the three-dimensional grids. The SRB grid was not axisym- 
metric because points were clustered in the small clearance between the ET and the 
SRB as shown in Figure 2.8c. 

The composite grid, consisting of the ET, SRB and ORB grids, is shown in 
Figure 2.7 at the plane of symmetry. The ET grid is the major grid and extends 
all the way to the far field boundary, while the ORB and the SRB grids are smaller 
computational domains with stings extending to the outflow boundary. Also visible 
in the figure is the hole boundary in the ET mesh cut out by the ORB. In Figure 2.8, 
the cross sectional view for all three component grids are presented at a constant £ 
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location. The hole boundary in each mesh, ET, SRB and ORB, is clearly visible and 
the hole points are removed due to the existence of the solid bodies. The values of 
the flow variables on the hole boundaries of the ET mesh are provided by the ORB 
and the SRB grids, while the flow variables on the hole boundaries of the ORB is 
updated by information from both the ET and the SRB grids, depending on the 
locations of the boundary points. The hole boundary of the SRB is entirely updated 
by the ET grid since the hole is cut out due to the presence of the ET only. The 
outer boundary update of the ORB grid, for the most part, comes from the ET grid, 
with a small portion near the SRB coming from the SRB grid. The outer boundary 
of the SRB grid is similarly updated — partly by the ET grid and partly by the ORB 
grid. Figure 2.9 shows the three-dimensional view of the hole cut out in the ET grid 
due to the ORB and the SRB. 




Figure 2.3: Detail views of the space shuttle vehicle 
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Figure 2.4: 


Simplified configuration and surface grid point distributions for prelim- 
inary supersonic flow calculations 



Figure 2.5: Improved configuration and surface grid point distributions 




Figure 2.6: Various computational planes of the orbiter grid 
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Figure 2.7: Symmetry planes of all grids 
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(a) ET grid 


(b) ORB grid 



(c) SRB grid 


Figure 2.8: Grid cross-section showing holes 
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Figure 2.9: Hole boundaries of ET grid 
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3. GOVERNING EQUATIONS AND NUMERICAL METHODS 

3.1 Introduction 

The thin-layer Navier-Stokes, Euler, and Poisson 1 equations were used to solve 
for the flowfields considered in this research. First used by Pulliam and Steger (1978), 
the thin-layer Navier-Stokes equations are based on the observation that for high 
Reynolds number flows, the viscous effects are confined to a thin-layer near rigid 
boundaries. The gradients in this layer vary very rapidly only in the direction nor- 
mal to the surface. Thus, all the viscous terms in the other two directions are dropped 
in this approximation. Coupled with the use of body-conforming grids and clustering 
of the grid lines near the body surface, the thin-layer Navier-Stokes equations can be 
used to properly resolve the flowfield around the body surface at high Reynolds num- 
bers. The Euler equations are derived from the same set of equations by dropping all 
the viscous terms. Finally, the Poisson equations are used to obtain a rough estimate 
of the flowfield. This solution is then fed back into the Navier-Stokes equation solver 
as the initial guess to improve the overall rate of convergence. It is believed that a 
great deal of CPU time is spent in damping out the initial transient for the flow solver. 
Thus, if a good initial guess is available, the overall convergence will be improved by 

iThe Poisson equation referred to throughout this report is a special case of the 
vorticity- velocity algorithm. 
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reducing the time in the transient state. In the following sections, the Navier-Stokes 
equations will first be presented followed by the thin-layer Navier-Stokes equations, 
and the vorticity-velocity algorithm. 


3.2 Navier-Stokes Equations 


The three-dimensional unsteady Navier-Stokes equations in the Cartesian coor- 
dinates can be given as (Peyret and Viviand 1975): 


dQ dE dF OG dEy dFy dGy 

dt dx ^ dy dz dx ^ dy dz 


where Q is the vector of the flow variables, E , F, and G represent the inviscid fluxes 
and Ey, F v and Gy correspond to the viscous fluxes. 



Q = 


pu 

pv 

pw 

e 


(3.2) 
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puv 


puw 
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pvw 
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u(e + p) 


V (e + p) 


w(e + p) 
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with 


T X x — A ( Ux + Vy + Wz'j + 2flUx 

(3.5) 

Tyy = A ( Ux + Vy + ®z) + 2/iVy 

(3.6) 

T ZZ = A ^ Ux + Vy + + 2flWz 

(3.7) 

Txy = r yx — M ( u y + v *) 

(3.8) 

Tjz = Tzx = M ( u z + w a;) 

(3.9) 

rt/z = r 2 y = (vz + wy) 

(3.10) 

fix = 7 kPt~^ dx^J + + VT X y + WTxz 

(3.11) 

fly — ~fnPr~^ dyej + uTyx + VT yy + wryz 

(3.12) 

fl z = 7 kPt~^ dz^l + ut Z x + VT zy + wtzz 

(3.13) 

ej = ep~ l - 0.5 (u 2 + v 2 4- w 2 ) 

(3-14) 


The Cartesian velocity components u, v, and w are nondimensionalized by the free 
stream speed of sound, aoo? whereas the density, p, and the total energy, e, are nondi- 
mensionalized by the free stream density, po ©? and poo^oo? respectively. Pressure 
be obtained from the perfect gas law: 


can 


p = ( 7 - 1) [e - 0.5 p (u 2 + v 2 + w 2 )] 


(3.15) 
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The ratio of the specific heat, 7 , is set equal to 1.4. k is the coefficient of thermal 
conductivity, p the dynamic viscosity, and A from the Stokes hypothesis is —2p/Z. 
The Reynolds number is Re and the Prandtl number is Pr. 


3. 2. 0.1 Generalized coordinates Body- fitted coordinates are employed in 
the numerical simulation to simplify the treatment of arbitrary geometries, espe- 
cially the imposition of boundary conditions. The flowfield is mapped onto a uni- 
formly spaced computational domain, and the transformed equations are maintained 
in strong conservation law form for the purpose of shock capturing. The generalized 
coordinate transformation is defined by 


T — t 

( = £ 

V = v{x,y,z,t) 

C = C 


(3.16) 


The transformed Navier-Stokes equations are given by: 


d T Q + d^(E- E v ) +d v (F- F v ) + d ^ (G -G v ) = 0 (3.17) 


p 


P U 

pu 


puU + £ x p 

pv 

, E = J -1 

pvU + (yP 

pw 


pwU + izP 

e 


(e + p)U - i t p 


(3.18) 
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P v 


pW 

puV + Tjxp 


puW + CxP 

pvV + Tjy p 

, G = J~ l 

pvW + (yp 

pwV + rizp 


pwW + C, z p 

(e + p)V - T ] t p 


(e + p) W - C t p 


and 

u = it + ( X u + ZyV + ( Z W 

V = 1)1 + 7 ] X U + TJyV + TJzW 
w = C t + Cxu + c yv + CzU) 

U, V, and W are contravariant velocity components. The viscous terms 


0 


E v = J~ l Re~ l 


ix^xx + iy T xy + £z r xz 
£x r yx + iy r yy + £z r yz 
ix^zx + £ y T zy + iz r zz 
ixflx + iyfiy + £zfiz 


0 


F v = J~ l Re~ l 


Vx^xx + ’HyTxy 4- fJzTxz 
Vx T yx + Vy T yy + Vz T yz 
Vx T zx 4- Vy T zy + Vz T zz 
Vxflx + VyPy 4- PzPz 


(3.19) 


(3.20) 


(3.21) 


(3.22) 
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G v = J~ 1 Re ~ 1 


(x T xx + iy T xy + (z T xz 
(xTyx + iy T yy *b Cz T yz | (3.23) 

CxTzx + Cy T zy + iz T zz 
ix0x + i yPy + Czfiz 

where the components of the shear-stress tensor are given in Eqs. (3.5-3.10). The 
respective Cartesian derivative terms are expanded according to the chain rule of 
partial differentiation; for example, 


u x — £x u £ + 7 x ^77 4- ( (3.24) 

The metric terms are obtained from the chain-rule expansion of x £, y^, z^, etc., and 
solved for iy, £z, etc., to give 


(x = j (yyZ^ - Zyy^j 

yx = J ( Z £V£ - 

iy - J ( z y x ( - x y z () 

Vy = J (x£Z£ - X( -z^ 

iz = J (x v y^ - Vtjx^) 

Vz = J (y£X£ - y( x {) 

Cx = j ( z v y £ - z^yTj) 

it = ~ x t£x - yriy ~ Z T ^ Z 

H 

1 

p- 

II 

rjt = -x t tjx - yrVy ~ z T rj z 

P- 

H 

Pi 

1 

SP 

Pi 

II 

it = ~ x tCx ~ yriy - z T Cz 


and the Jacobian of the coordinate transformation is given as 

J = (*{VtiZ£ + x^y^z v + Xjjy^ - x^y^zy - x v y^z^ - x^z^J _1 (3.26) 


The metric terms, iy , etc., are all differenced using central differencing for the 
interior points, and second-order one-sided differencing for the boundary points. How- 
ever, in rare cases when the spacing between ( = 2 and ( = 1 is much smaller than 
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that between < = 3 and C = 2, the second-order one-sided differencing in the < deriva- 
tives causes the Jacobian to become negative. In these cases, switching to first-order 
one-sided differencing remedies this problem. 

3.3 Thin-Layer Navier-Stokes Equations 

The numerical procedure used for the thin-layer Navier-Stokes equations in this 
research was developed by Steger et al. (1986) in a program called F3D. For high 
Reynolds number flows, the viscous terms, E v and F v , can generally be neglected 
based on the same arguments used for the boundary-layer approximation. The cross 
derivative terms in G v are also dropped for the same reasons. However, in this 
research, the viscous term, F v , was treated the same way as the viscous term, G v , so 
that viscous effects could be accounted for in either the £ or q directions or both. The 
remaining of the viscous terms, F v and G v , were collected into the right-hand-side, 
i.e., Re- 1 (d^S + d v R), of the thin-layer Navier-Stokes equations as shown below. 

drQ + d^E + d v F + d^G = Re _1 ( d v R + d^S) (3.27) 

0 

/i + Vy + Vz) U V + (/*/ 3 ) (v^ u V + Vy v V + VzWy) Vx 

H (i]l + Vy+ Vz) v V + (m/ 3 ) {vx u v + Vy v q + VzWtj) Vy 
R = J - 1 fl(vx + Vy + Vz) W V + W 3 ) (v x U V + Vy v V + VzWy) Vz 

{(vl + Vy + V%) 

x |o.5/i (n 2 + v 2 + w 2 )^ + fiPr- 1 (7 - 1) _1 )v] 

+ (/i/3) (rjxu -1- r ]yv + T] Z w) + Vy v V + Vz^y) } 
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0 

M (Cx + Cy + Cz) u ( + (m/3) (Cx«£ + Cy*>£ + Czio() Cx 
M (Cx + Cy + Cz) v £ + (m/ 3 ) (Cx«£ + Cy v ( + Cz w () Cy 

5 = /» (cl + cj + Cz 2 ) ^ + W 3 ) (<* u c + ^c + CzW c) Cz (3 - 29) 

{ (Cx + Cy + Cz) 

x ^0.5/r 4- w 2 )^ + (7 — 1)~1 (« 2 )^] 

+ (/i/3) (Cxu + Cyv + Czw) (Cx«£ + CyV( + Cz«^) } 

In the implementation of the thin-layer Navier-Stokes equations, the following equa- 
tion was used instead of Eq. (3.27). 

dr (Q - Qoo) + d^E - Eoo) + d-q {P - F 0 o) + d^(G - Goo ) ^ 

= Re- 1 (d v R + d^S) 

The “free stream subtraction” used in the above equation is to avoid the errors 
introduced from the approximations in computing the metric terms x^, y^, etc. The 
Roo and §00 are neglected since they are small for high Reynolds number flows. An 
implicit two-factor approximate factorization scheme, shown below, that uses central 
differencing in the 77 and ( directions and upwind differencing in the £ direction is 
used for the above equation. 

I + h8 b {A+) n + h8^C n - hRe~^8 < *J—^M n J - 0^] 

x\l + hs{(A-) n + h8 v B n - hRe~ 1 8 v J- l N n J - D : |J A Q n = 

C v ' 1 '1 (3.31) 

-Af{^(£+) n + d{E~) n + 5 v F n + ^<$ n - Re~^8jjR n -Re^S^S 71 } 

—{Delri + D e \()Q n 

where 8 b and 8f are the backward and forward difference operators, 8 the central- 
difference operator, 8 the mid-point central difference operator for the viscous terms 
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and .Dj and D e the implicit and explicit smoothing operators used in the centrally dif- 
ferenced directions, V and (. The scheme is capable of achieving first-order (h = At) 
or second-order time accuracy (h = At/2). The matrices A ± , B, C, M, and N are 
the results from the local linearization of the fluxes, £ ± , F, G, R, and S about the 
previous time level, e.g., 

jP n +l = F n + B n AQ n + 0 (h 2 ) (3.32) 

The implicit approximately factored scheme, Eq. (3.31), is solved by marching in the 
£ direction using two sequential sweeps of the block tridiagonal inversion procedure, 
one in the t / and the other in the ( direction, at each constant £ plane. 

In the chimera scheme, the flow solver developed for a single curvilinear grid has 
to be modified to account for the hole points (including the hole boundary points) 
introduced in the overset grid. At such points, the values should be kept unchanged 
since these points are either at the hole boundary or within a body (or a user specified 
boundary zone), and are updated by other grids or assume no meaningful values. An 
array of values i b , 1 at the regular points, and 0 at the hole points, is thus introduced 
into the approximately factored scheme to turn off the finite differencing at the hole 
or hole boundary points. The following shows how i b is used in the differencing 

scheme. 

[/ + i b (h8 b (A+) n + h8 c C n - hRe- l ~8 c J- 1 M n J - I?;I C )] 
x \l + i b (h^(A~) n + h8 v B n - hRe~ 1 8 r] J~ 1 N n J - = 

-i b [A t8 b (E+) n + 8f(E~) n + 8-q F n + 8^ 

-Re- l 8 v R n - Re^^S 71 4- (D e \rj + Del^Q 71 ] 


(3.33) 


40 


Thus, for hole points, i ^ = 0, and the above equation reduces to 


AQ n = 0 


(3.34) 


and Q remains constant. If only three-point central differencing is used, the scheme 
requires no further modifications. However, when differencing the points adjacent 
to a hole boundary, finite-difference operators that require information beyond the 
adjacent points on either side of the differenced point will need to be modified since 
only the hole boundary points are updated from other grids, and the hole points 
do not contain meaningful data. Therefore, difference operators adjacent to a hole 
boundary will need to be modified such that data from the hole points are excluded 
from the calculation. The same i ^ mentioned above can be used to achieve this 
purpose. For example, the right hand side dissipation term in the £ direction is 
currently implemented as 


- (S±4±a) W /+1 - - «/) 

7 VAQ/-VAQ f _ 1 \ _ / qf+a/.n ' 

y l+aj+aj_-^ J y 2 J 'Ql l) 

with e = 0(.l), VA a second-order differencing, and 

a = / 1 + m qq x ( n± i ~ 2 jn +pi-i \ 

16 \P/+1 + 2pi +Pi-i) 


Here g is a modified spectral radius of the matrix C (see Eq. 3.31), 


(3.35) 


(3.36) 


Q — |Ca: u + Cy v + (zH + \j (Cx^ + Ct/^ + Cz^)c^ + .01 


(3.37) 
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2 2 2 \ 

with c 2 = 2£(1 _ 0) + + w _i and 0 < /? < 1 where the choice fi = 1 

reduces the smoothing in the boundary layer. The parameter a determined from the 
pressure gradient is used to switch from second-order to fourth-order smoothing. To 
avoid using data from the hole points, the differencing is modified as 

= (5) ($jtl + $0 x 

- (^) (o« + i - «w] <3 38) 

-(2) + 7731) x 

(Ql - <?,_,) 

Thus the fourth-order differencing is reduced to an uncentered second-order differ- 
encing adjacent to a hole boundary point. 

3. 3. 0. 2 Boundary conditions Explicit boundary conditions were used for 
the thin-layer Navier-Stokes equations. Each boundary condition was coded in a 
separate subroutine; thus, additional boundary conditions can be conveniently added 
as need arises. The boundary conditions implemented include inviscid/viscous wall 
conditions, far-field conditions, axis conditions, wing cut conditions, symmetry plane 
conditions, periodic conditions, and overset grid hole and outer boundary conditions. 

3. 3. 0. 3 Wall conditions For viscous flows, the no-slip boundary condition 

is enforced by setting the velocities on the wall boundary to be zero. For inviscid 
flows, the tangency boundary condition is implemented by setting the contravariant 
velocities = 0; in other words, the fluid flow is not allowed to go through the 
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wall boundary. The values of the contravariant velocities U and V are extrapolated 
from the interior points at ( = 2 and ( — 3. The Cartesian component of the velocities 
can be found by solving Equation ( 3.20). The density on the wall is obtained through 
zeroth-order extrapolation from the interior points at ( = 2, i.e., p^~i ~ P(=2' This 
is a reasonable approximation in the current study since the grids used are usually 
clustered near the surface and the grid spacing there is very small. 

The pressure on the surface is obtained from a normal momentum relation found 
by combining the three transformed momentum equations (Pulliam and Steger 1978), 

(&cCx + (y(y + (z(z) P£ + ( Vx(x + Vy(y + VzCz^ Prj 
+ (Cx + Cy + Cz) P( ~ (3.39) 

—pU ((xu^ + (yv^ + (zw^) — pV (^Cx^tj + (y v rf + (zWrf) 

The above equation is solved using central differencing for the ( and rj derivatives, 
second-order one-sided differencing for the ( derivatives, and values from the interior 
points at ( = 2 and ( = 3. For viscous flows, U = V = 0 is used in the above 
equation. 


3. 3.0.4 Far-field conditions At the far-field boundary, two different bound- 
ary conditions were implemented. The first one is the Dirichlet boundary condition 
which simply sets the free stream values for points on the ( = (max surface. Gener- 
ally, this boundary condition requires the largest computational domain among the 
possible far-field boundary conditions, and is valid only when the far-field boundary 
is far enough from the solid boundary. 


Q far-field ~ Q°° 


(3.40) 




Figure 3.1: Surface grid for an ellipsoid 


The other far-field boundary condition implemented is the outflow boundary 
condition. It assumes zero gradient at the downstream boundary. Of those grids 
implemented with this boundary condition are the ET, SRB and ORB grids. 

Q f =Q f , (3-41) 

^ Kmax ^Kmax-l 

3. 3. 0.5 Axis conditions For the grid extending from both ends of an ellip- 
soid (Figure 3.1) or from the ET nose, the grid points fall onto a single line which 
is referred to as an axis for convenience. The values of the variables on the axis are 
evaluated from the average of the variables at the neighboring points surrounding the 
axis. This condition was implemented for £ = constant and £ = constant axes. 

««i. = l t V < 3 - 42 > 

1 = 1 

Note that this averaging should not include hole points since they contain no mean- 
ingful data. Other geometries, like the SRB, ORB, ellipsoid, and sphere, also used 
this boundary condition at the axis. 
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3. 3. 0. 0 Wing cut conditions The orbiter grid was generated with a wing 
sting extending from the trailing edge to the outflow boundary. For the flow simula- 
tion to more accurately mimic the real flowfield, a condition was introduced to allow 
the flow to go through the wing sting without regenerating a new sting-less orbiter 
grid. What this condition did was to set the flow variables at the top and bottom 
surfaces (at the same span wise location) to equal values. Since the thickness of the 
sting is very small, this is considered a reasonable approximation. 

3. 3. 0. 7 Symmetry plane conditions To save computer time, a symmetry 
plane condition was used as the situation permits, as with the ET and ORB calcu- 
lations. This condition assumes no penetration of fluid flow through the symmetry 
plane and sets values for other variables at points on opposite sides of the symmetry 
plane equal. 

3.3.0. 8 Periodic conditions This condition ensures that flow variables at 
points on the planes at 7] = 1 (0°) and j] — 7]max (360°) are of equal values. This 
condition was used for the SRB grid, where the flow had to be solved for a complete 
360° circle. 

Qr]—l ~ Qr\~max (3.43) 

3. 3. 0. 9 Overset grid hole and outer boundary conditions As was ex- 
plained in the chapter on “Chimera Approach and Grid Topology”, the interpolation 
is required to transfer information (through the hole and outer boundaries) among 
the different grids of the composite geometry. For this study, trilinear interpolation 
was used. 
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3.3.1 Solution Procedure 

As described in the chapter on “Chimera Approach and Grid Topology”, the 
major advantage of the chimera approach is the embedded flexibility in the chimera 
scheme which allows one to blend different grids, different flow solvers, etc. Here, 
a flow solver developed for a single curvilinear grid can easily be tailored for the 
chimera scheme by just adding a control loop outside the major flow simulation loop 
to update the grid interface boundaries and provide ways to blank out the hole and 
hole boundary points. 

In the current implementation, data from each grid and the boundary interface 
arrays are brought in sequentially from external storage (high speed disk or CRAY 
solid state device, SSD) at each time step. The flow variables on the hole and outer^ 
boundaries are first updated by the values stored in the boundary interface arrays. 
Then the solution is updated by the flow solver and the imposed boundary conditions. 
The boundary interface data that the current grid sends to other grids are then loaded 
into the boundary interface arrays and all the arrays are sent back to the external 
storage. The next grid is then brought in and so on. 

3.4 Vorticity- Velocity Formulation 

The vorticity- velocity formulation consists of: 

1. Poisson equations — derived from the continuity equation. 

2. Crocco relations — derived from the momentum equations for inviscid 
flows and the Gibbs function, a thermodynamic relation based on the first 

^Outer boundary is updated only if the current grid is embedded in other grids. 
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and second laws of thermodynamics. 

3. Convection of entropy and stagnation enthalpy 

4. Bernoulli equation — derived from Crocco relations and the equation of 
state. 

5. Vorticity consistency condition — a vector identity which enforces the 
conservation of vorticity in the flowfield. 


The details of each equation will be described in the following sections in the order 
listed above. 

The vorticity- velocity formulation is a set of weakly coupled equations. Hence, 
it can be solved more efficiently than the Navier-Stokes or Euler equations since only 
scalar tridiagonal (rather than block tridiagonal) systems of equations are involved in 
the computation. However, this algorithm won’t be as accurate as the traditional flow 
solver using Navier-Stokes equations because viscous terms are omitted. Currently, 
the algorithm is limited to subsonic flows and is intended to be used as a diagnostic 
tool for testing grids (especially, in the case of complex geometry consisting of overset 
grids) or to offer an approximate solution as an initial guess for the more accurate 
but more expensive solvers using the Navier-Stokes equations. In the future, the 
current algorithm may be extended to treat viscous flows by combining with the 
boundary-layer or Navier-Stokes equations. 
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3.4.1 Poisson equations 

The Poisson equations are derived from the continuity equation. First, by writing 
the continuity equation in non-conservative form, 

V- V + 1 ? = 0 (3.44) 


where 


* = ?£■■ v 


(3.45) 


and then via differentiation of Eq. (3.44) and the use of the following vector identity, 


V x (Vx V) = V(V- V)-V z V 


(3.46) 


the vector form of the Poisson equations are obtained. 


V 2 V +V x (Vx V) + Vtf = 0 


(3.47) 


where Vx V can be replaced by the vorticity vector, u ;, 


uj = ( Wy — V Z )i + ( Uz — U>x)j + (v X ~ «t/)£ 


(3.48) 


and thus 

Vx(VxV) = V x u> 

= ( w 3 y ~ “2 z)» + ( u lz ~ “Zx)] + ( w 2x “ w ly)* 
For three-dimensional flows, the Poisson equations can be written as 

V 2 u + tix + <^3 y - “2 z = 0 
V 2 u + tiy + u\ z - U> 3a . = 0 
V 2 W + + U>2 X -H>ly = 0 


(3.49) 


(3.50) 

(3.51) 


(3.52) 
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where is the Laplacian operator, = dxx +dyy + dzz • 

For irrotational flow, the 

vorticity terms can be dropped from the Poisson equations, 

resulting in the following 

equations: 


U XX + Uyy + Uzz + — 0 

(3.53) 

Vxx 4* Vyy 4* Vzz 4- $y = 0 

(3.54) 

Wxx 4* w yy 4 _ tt Jzz "1" $ z = 0 

(3.55) 


3. 4. 1.1 Generalized coordinates The Poisson equations in the generalized 
coordinate system are derived using the chain rule of partial differentiation and the 
vector identities resulting from the coordinate transformation. Only the results are 
shown here; for details of the derivation, see Appendix A. 

+ 

V , (3-56) 

+ 4 = 0 

c J 




fu — txdg' + +£y w 3 £ + r ly u) Zri +.Cy u; 3£ — — Vz<^2t] ~ Cz^C (3-59) 

fv = + + £z c * ; l£ + Vz<^Itj + Cz w i£ ~Cx^( (3-60) 
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fw = ud(+ r izdTi+(zd ( ;+ZxU2z+ r ixu2 V +( x U2c~ty UJ H~ 11 y ul iv~ < 'y u} K ^ 3 ' 61 ) 


and 

a ££ = (x£x + £t/£y + iziz (3.62) 

<27777 = VxVx + 7 ly T \y + VzVz (3.63) 

= CxCx + CyCy + CzCz (3.64) 

a ^ = ixVx + (yVy + (z^z (3.65) 

cl££ = £xCx + £yCy + (z(z (3.66) 

= t/xCx + VyCy + ^zCz (3.67) 


3. 4. 1.2 Finite-difference formulations The cross derivative terms {Q^£, 

Qtf, Qfr’ Qcr <?£<> and ®v0 are 311 lagged at the previous time level; thus the 

finite-difference formulation of the Poisson equations can be written as 

+ 8 ( (”^ 6 c)] ^ n+1 " 2 

= - + -$^<) + 8i i {~% L8 t + ~ J f 8 () (3 ’ 68) 

where <5 designates the central difference operator and h the time step size. 

Note that an additional time derivative term, (Q n+1 -Q n )/h, is added to the 
Poisson equations to allow the use of an approximately factored scheme, which will be 
explained later. This added time derivative term should not affect the solution when 
convergence is reached, since Q n — Q n+1 at convergence. The above equation can 
be rewritten in delta-form (treating the change (<? n+1 - Q n ) as the unknown vector) 
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as shown below. 


{' - h K (^ f {) + b (Tb) + b (^* c )]} («" +I - <?") 

= hu, ( a K s ( +a (f ' +a v; s < 'j + s v 

<? n + hto£ 


+ S C ^ a « S < +a ^ +a vC Sl l 
= RHS 


(3.69) 


The relaxation parameter (time step), h, varies with the Jacobian of the grid point as 
shown in the equation below; since the Jacobian varies from point to point, values of 
h also vary at different locations. This procedure is called local time stepping since 
the time step size is determined by the local grid size. However, this method is only 
applicable if a steady state solution is desired as in the present study. 


, AT(1 + 0.005 J) 

fl = 

1 ~f J 


(3.70) 


A geometric sequence was also used to determine the relaxation parameter, h. 


h = 


V 

M. 


i—i ■ 

7V^I 


-1 


i = 1,2,3 ...JV 


(3.71) 


A second overrelaxation parameter, a>, for the source term, y, was used to speed up 
the rate of convergence. Values of Aj = 0.02 to 1, A 2 = 0.005 to 0.01, N = 4, and 
u> = 1.0 and 1.8 have been used in the course of the study. If Aj and A 2 are chosen 
properly, the geometric sequence method for determining the relaxation parameter, 
h , can work just as efficiently as the local time stepping. 

The three-step approximately factored scheme, as shown below, was used to 
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solve the Poisson equations. The scheme is second-order accurate in space. 


■cW 


1-hS 


[l - h6 v (“J 2 ) s v] x 

(<5 n+1 - Q n ) = RHS + 0{h 2 ) 


(3.72) 


Central differencing is used for all the terms in the above equation. For chimera 

O 

overset grids, the hole and hole boundary points should be left unchanged 0 ; thus, the 
blanking array, needs to be introduced into the above equation to shut off the 
differencing scheme for these points. Simply by replacing h with i^h, and RHS with 
ifoRHS, the Poisson equations become 


( Q n+1 -Q n ) = 0 


(3.73) 


for hole or hole boundary points ( i ^ = 0). Thus, values of flow variables for these 
points remain unchanged. 

The following shows the finite- difference formulation used for the Poisson equa- 
tions. Unless stated otherwise, indices, j, k , and l, used in finite-difference expressions 
refer to grid points in the £, 77 and ( directions, respectively, and throughout this re- 
port, a subscript in a finite-difference expression is not shown unless it varies, e.g., 

Qj + 1 = 

( = ( Q i+}~Qj) 

v J h ^ (3.74) 

2 Sf 

^Hole points are points which lie within a body or user specified boundary zone 
and should not be solved, while hole boundary points, considered as part of the 
boundary, are updated by other grids and should not be solved either. 

^See Section 3.3 “Thin-Layer Navier-Stokes Equations” for details. 
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The finite- difference expressions for ^ ^ and ^ ^ ^ can be obtained by 

replacing the j index above with the corresponding k and l indices. For the cross 
derivative terms, the finite-difference expressions are given in the form: 


W, ■* |(^L - (*£)„ 


1 




(Qj+i 1 fc+i~Qi+i,fe-i) 




( 3 . 75 ) 


(^) ^ k+l-Qj-lJt-l) 


Again by replacing the above j and k indices with the corresponding indices for 
other cross derivative terms, p£i) , (^) , (=s£t) , 

( a„sQv \ 

and 


c 


, the respective finite-difference expressions can be obtained. For 
the chimera overset grids, the blanking array, needs to be incorporated into the 
finite-difference expressions for points involving the flow variables at the hole points. 
Only the cross derivatives need to be modified since only those terms involve flow 


variables at hole points. Given below is the modified finite-difference expression for 

'au 

T , 




(*&) - 


( %+!,&+! +* 6 j + 1 ,k-l +t bj - l,fc + 1 l,fc + 1 


a 


* 6 j+l,fc+ll "7 



J +l bj+l,k-l 


UtiMh ' 




( 3 . 76 ) 


+i bj-l,k+l 


a 








53 


Central differencing was used for all the Q derivatives in the above equation. For other 
cross derivative terms, the finite-difference expressions can be similarly obtained. The 
resulting finite-difference expressions for the Poisson equations can be given in the 

form: 

L^LrfL^{Q n+1 ~ Q n ) = R ( 3 - 77 ) 


The solution procedure for the above equation involves three sequential ADI sweeps 
in the (, i] and < directions. The following illustrates the step-by-step procedure for 
this algorithm. 


L^AQ* 

-- R 

LrjAQ** 

= A Q* 

L^AQ n+l 

= A Q** 

QU+1 

= Q n + AQ 


n+1 


(3.78) 


3.4.2 Crocco relations 

For inviscid flows, the momentum equation can be written as 


p(V 



(3.79) 


Substituting the vector identity below 

(V -V) V= V 



V xw 


(3.80) 


into equation(3.79), the inviscid momentum equation can then be rewritten as 


V 



— ► _ Vp 

- V xu = - 

P 


(3.81) 


Crocco’s equation is then obtained from the above equation by expressing the pressure 
gradient in terms of an entropy and an enthalpy gradients using the Gibbs equation 
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(3.82), which is a thermodynamic relation derived from the first and second laws of 
thermodynamics . 

TVs = Vh~ — (3.82) 

P 

For a steady, adiabatic flow (H = h + = constant), the Crocco relation can thus 

be written as 

V xu = -TVs (3.83) 

3.4.3 Convection of entropy and stagnation enthalpy 

By taking the dot product of the velocity vector and the Crocco equation (3.83), 
and assuming constant total enthalpy, the convection of entropy is found to be 

V -Vs = 0 (3.84) 

or in Cartesian coordinates 

usx + vsy 4- wsz — 0 (3.85) 

By applying the chain rule of differentiation to the above equation, and the defi- 
nition of contravariant velocities in equations (3.20), it can be shown that in the 
computational domain, the entropy convection can be expressed as 

Use + Vs v + Ws^ = 0 (3.86) 

Similar equations can be obtained for convection of stagnation enthalpy by substi- 
tuting stagnation enthalpy for entropy. 


55 


3. 4. 3.1 Finite-difference formulations To illustrate the differencing used 
with the convection equations, consider the entropy convection equation in gener- 
alized coordinates, equation (3.86). Three-point second-order central differencing is 
used in 77 and £ and upwind differencing is used in £. Thus U s g is differenced as 


Use = 


u+m a. , v-w\j 


Sfj + 


S J £ s 


(3.87) 


Depending on the values of the contravariant velocity, U, either forward or back- 
ward upwind differencing is used. The convection equation is then solved using the 
approximate factorization algorithm used by Bridgeman et al. (1982). 

(I + h U+ ^ S h £ + hW8(-){I + h U ~^ U h ^ + hV6rj)(s n+1 - s n ) = 

(3.88 ) 

+ V*n + W8(-)s n 

where h is a relaxation parameter, h > 0. Adding second-order numerical dissipation 
in the 77 and £ directions to the above equation gives 
(/ + hU + 6% + hW8 ^ - h|W| AV|^)x 

(I + hU~8 ? + hVSrj - /i|V|AV| T? )(3 re + 1 - s n ) = (3.89) 

-h [u+6% + US £ 4- V 8rj + W8 £ + |W'|AV|, 7 + |V|AV| C ] s n 


where u+ m and v - m Hzfl. 


3.4.4 Bernoulli equation 

By assuming that the flow is steady, inviscid, adiabatic, and using the perfect 
gas relations, the Bernoulli equation can be derived from the Crocco relations and 
perfect gas relations (Anderson et al., 1984): 

1 

e -(s-soo)/R (3.90) 


Poo 


1 + 


7 - 1 


2 , 2 , 2 ' 

„,2 u + v + w 

Moo 2 

a oo 


7- 


56 


For irrotational flow, the entropy correction term, e 9 oo)/R ? can b e dropped. 

To avoid expensive exponential operations, the exponential term was expanded 
using the binomial expansion: 


n 


(1 ± a) = 1 ± an H — a: ± 


(n-1) 2 , n(n- l)(n - 2) 3 


2x3 


oc -f- . . . (3.91 ) 


where 


Written for chaining 


7 _ 1 / „ r 2 

a = — — I o 


a oo 


(3.92) 


(1 + a) n = 1 + an(l + a — - ^ (1 + a ^ n - 2 ^ (1 + ^ . . .))) (3.93) 


or 


or 


(1 + a) n - 

1 + a(n + + a n(n-l)(n-2)(n-3) })) 

(1 + a) n = 1 + «(ci + «(c2 + 01(03 + <*04))) 


where 


C 1 


7 - 1 


^ = *> 
=3 = f(^r-2) 
c 4 = ^^r- 3 ' 


(3.94) 


(3.95) 

(3.96) 

(3.97) 

(3.98) 


(3.99) 



57 


3.4.5 Vorticity consistency condition 
From the vector identity 

V • (Vx V) = V • w = 0 (3.100) 

the consistency relation is obtained 

+ dy <*>2 + 3 = ® ( 3 . 101 ) 

This consistency condition is used with the tangency boundary condition to form the 
boundary condition imposed for the solid boundary. 

3.4.0 Boundary conditions 

Since the boundary conditions coded for the thin-layer Navier- Stokes equations 
are of the modular form, they are readily available for the vorticity- velocity formula- 
tion. The only boundary condition implemented specifically for the vorticity- velocity 
formulation is the tangency/vorticity consistency condition. This boundary condition 
was implemented implicitly to overcome the slow convergence found in the course of 
this research. 


3. 4. 6.1 Tangency/vorticity consistency boundary condition The tan- 
gency (or no-flow-through) condition was imposed on a solid boundary by setting the 
contravariant velocity to zero in the direction normal to the solid boundary. This 
was combined with the definition of vorticity to form the boundary condition. The 
following shows the derivation of this boundary condition. From the vorticity defini- 
tion 

(yW£ - (zV( = f\ = - {tyw^ + TJyWy - U v ( ~ VzVjj) (3.102) 
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Cz u £ — Cx^£ — $2 — ^2 _ (Cz u £ + t Jz u t] — £x w £ ~ t )x w t]) (3.103) 

Cx^£ - C y u £ = fz = ^3 ~ (Cz v C + VxVrj — (yu^ — Tj y u v ) (3.104) 

differenced on the wall as 

(yw - Cyw* - ( z v + C zv * = -AC f± (3.105) 

Cz« - Czu* - (xu> + Cxw* = -AC /2 (3.106) 

Cxv - Cxv* - (yu + Cy u * = -AC /3 (3.107) 


where u*,v* and w* are values of u,v, and w at A C above the wall. The tangency 
condition on the C = 0 surface is given by 


Cx u 4 " (yv + Czw — 0 


(3.108) 


Rewrite the above equations in matrix form: 
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(3.109) 

Solving these four equations for the three unknowns'* u, v , and w on the surface by a 
generalized inverse (or using the vorticity relations to remove the other components 


^The first three equations in Eq. (3.109) defined by the vorticity definition are not 
linearly independent; therefore, Eq. (3.109), in essence, consists of three equations 
and three unknowns. 
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from the tangency relation, or by multiplying the above equation by the transpose 
of matrix A) gives the vorticity/tangency relations: 
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(3.110) 


where l 2 = Cx + Cy + Cz 

The tangency boundary condition was solved implicitly in the ( direction along 
with the Poisson equations to overcome the tendency for slow convergence due to the 
use of an extremely fine grid spacing near the solid boundary. The details of this 
procedure are described in Appendix B. 


3.4.7 Solution procedure 

3.4.7. 1 Irrotational flow Since the Poisson equations are weakly coupled, 
the Thomas algorithm can be used to invert each tridiagonal system of equations 
sequentially without resorting to the use of the more expensive block tridiagonal 
solver. At each grid point, an initial guess for u,v,w, and p is made. The following 
procedure is then implemented to obtain the flow solution. 

1. Velocities u,v,w Sequentially update u, v, w from the Poisson equa- 


tions. 
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2. Density, p Update density from Bernoulli equation using the most 
recently calculated values of u,v and w. 

The above procedure is repeated until convergence. The solution was considered 
converged when the L2-norm dropped three orders of magnitude. 

3. 4. 7. 2 Rotational flow At each grid point, an initial guess for u, v, w, p, 
s, H, u>i , u> 2 , and u>g is made. The following procedure is then implemented to obtain 
the flow solution. 

1. Velocities u,v,w Sequentially update u,v,w from the Poisson equations 

for assumed values of and 0/3. 

2. Entropy, s, and stagnation enthalpy, H Update s and H from the 
corresponding convection equations using updated values of u,v and w. 
Currently, H is assumed a constant. 

3. Density, p Update density from the Bernoulli equation using the most 
recently calculated values of u,v,w and s. 

4. Vorticities, u>i,u>2, and 0/3 Vorticity components are calculated using 
the Crocco relations and the consistency condition for vorticity. 


This procedure is repeated until convergence is reached. 
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4. RESULTS AND DISCUSSION 


The thin-layer Navier-Stokes solutions for the integrated space shuttle vehicle 
are presented first in this chapter followed by the solutions from the vorticity- velocity 
formulation for several geometries. 

4.1 Thin-Layer Navier-Stokes Solutions 

The grids used in the flow simulation for the integrated space shuttle vehicle were 
generated using a hyperbolic grid generator (Steger and Rizk, 1985). The ET grid 
extends all the way to the far-field boundary where free stream values were assumed 
while the SRB and ORB grids only fill a much smaller portion of the computational 
domain so that duplicated effort in numerical calculations can be avoided. The holes 
cut out in each grid as seen in Figure 2.8 are usually larger than the dimensions of 
the embedded body and preferably larger than the dimensions of the boundary layer 
so that the chimera interpolation won’t be carried out in regions of high gradients. 
However, it may not always be possible to carry out the interpolation outside the 
boundary layer if two grids are very close to each other, like inside the clearance 
between the ET/SRB and ET/ORB. Thus, more points are needed in these regions 
to make the interpolation accurate or possible. The computational domain extended 
to the outflow boundary where zero gradients of the flow variables were implemented 
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for all three component grids (ET, SRB and ORB). In the circumferential direction, 
symmetry plane boundary conditions were used for the ET, ORB and the forward 
and aft attach grids while periodic boundary conditions were employed for the SRB 
grid. The no-slip wall boundary conditions were applied on the body surface for all 
grids. The axis boundary condition was implemented for axes extending upstream 
from the nose of the respective grid to the far-field boundary. Initially, a smaller grid 
of 289, 212 points was used in the calculation for a flow of free stream Mach number 
2.0, and gradually the grid size was increased to about 750,000 points including 
treatment of the idealized attach hardware between the ET and ORB for a transonic 
flow of Mach number 1.05. Other than the grid points used for the additional grids 
(forward and aft attach grids), the added grid points were concentrated in regions 
where the important flow physics was expected to occur, e.g., regions around shock 
waves, the recirculation region at the back end of the ET, and the regions near 
the attach hardware where the flow separates due to the obstruction of the attach 
hardware. 

Before proceeding to compute the flowfield around the integrated space shuttle 
vehicle, the Pegasus code (Benek et al., 1985) was used to obtain the chimera inter- 
polation data between the component grids. Then the process of obtaining the flow 
solution started from a uniform free stream flow with the wall velocities gradually 
reduced to zero in 30 iterations to minimize the effect of possible oscillations resulting 
from setting the wall velocities to zero too abruptly. The flow on the ET grid was 
computed first with the necessary boundary conditions — far-field, wall, outflow, and 
chimera boundary conditions (the flow variables on the hole fringe assume free stream 
values initially). The solution on the ET grid was then used to update the outer and 
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the hole boundaries of the SRB and ORB grids using the interpolation information 
provided by Pegasus. Then the solutions on the ET grid were written to a high speed 
solid state device (SSD) or on external disk. The solution on the SRB grid was then 
computed using the updated outer and hole boundary values and other necessary 
boundary conditions. As with the ET grid, the solution on the SRB grid was used to 
update the ORB outer and hole boundaries and the ET hole boundaries (no need to 
update the ET outer boundary since it was outside the SRB grid and assumed free 
stream values) before sending it back to SSD or external disks. A similar process was 
carried out for the ORB grid and other grids, if there were any. This completed an 
iteration for the entire overset grid. The next iteration was then repeated beginning 
with the ET grid and so forth. 

Flows of different free stream Mach numbers — 0.6, 0.9, 1.05, 1.55 and 2.0 
— were calculated. The Mach number at which the maximum pressure loading on 
the space shuttle occurs during ascent is usually close to this range. The Reynolds 
numbers used in the computation were taken from the wind tunnel tests, and the 
angles of attack from the actual flight data. To assess the feasibility of the chimera 
approach for the integrated space shuttle vehicle configuration, a supersonic flow 
testcase, Moo = 2.0, was chosen first since it is cheaper and easier to compute than 
the more difficult subsonic flow case. A coarse grid of 289,212 points was used for 
the entire overset grid. Many simplifications were made for this grid as can be seen 
in Figure 2.4. The stings behind the SRB and ORB were used to simulate the plume 
effects while the ORB vertical tails, the ET/SRB attach ring, and the ET/ORB 
attach hardware were all missing from this grid. Since the supersonic flow has only 
limited upstream influence, even with these simplifications, the numerical solution 
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was expected to capture some important flow phenomena. In Figure 4.1 the simulated 
oil flow obtained from PL0T3D (Buning and Steger, 1985) for a constant £ plane is 
compared with the oil flow from the wind tunnel test for the integrated space shuttle 
vehicle at Moo = 2.0 and a = —4°. The Reynolds number was not correctly modeled 
at the time of this simulation due to the lack of wind tunnel test data. Nonetheless, 
the simulated oil flow patterns show reasonable agreement with the wind tunnel 
oil flow results for portions of the geometry. The differences can be attributed to 
inaccurate modeling of the geometry, especially in regions near the ET/ORB attach 
hardware, the protuberances on the SRB surface, and the back end of the ET where 
the wind tunnel model does not have a sting. A quantitative comparison of the 
surface pressure coefficients for the orbiter is presented in Figure 4.3a. The flight 
data (Rockwell International, 1983) in the figure were taken on a constant angle 
station (</> = 70°) along the side of the orbiter fuselage. As noted in the figure, 
there are discrepancies in the angle of attack and the elevon deflection; however, the 
computation does show a trend similar to the flight data. 

For the Moo = 1-55 case, the same grid was used but with the correct wind tunnel 
Reynolds number, Re = 3.2 x 10® /ft, and the angle of attack, a = —6°. Figure 4.2 
shows the surface pressure coefficient comparisons between the computation and the 
wind tunnel tests (Spangler, 1981). The 3% wind tunnel model of the integrated 
space shuttle vehicle was equipped with 1538 pressure taps which provided enough 
data to allow extraction of meaningful pressure contours from the experimental data. 
The wind tunnel test was carried out by Rockwell International, Inc. For convenience 
of reproduction, gray scale contours, instead of color contours, of the surface pressure 
coefficient are given in Figure 4.2. Unlike the color contours, the gray scale contours 
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Figure 4.1: 


Wind tunnel and simulated surface oil flow for the integrated vehicle at 
Moo = 2.0 and a = —4° 
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of the surface pressure coefficient in Figure 4.2 do not show the magnitude of the 
pressure coefficient. However, the variational change does offer insight on the extent 
of the expansion and compression regions. Overall, the computation does agree 
reasonably well with the experimental results except where the computational model 
failed to model the real geometry, e.g., the lack of the ET/ORB attach hardware and 
the ET/SRB attach ring and the redundant ET sting. In Figure 4.3b, the comparison 
of surface pressure coefficients for the orbiter along the <f> = 70° line shows that the 
numerical solution is generally in agreement with the experimental and flight data. 
The major discrepancy is at the trailing edge of the orbiter wing where the numerical 
model did not correctly model the elevon deflection. 

For the flows of Moo = 0.6 and 0.9, a grid of 400,902 points was used. This 
grid still contains only the three major components of the integrated space shuttle 
vehicle, i.e., the ET, SRB and ORB. The added points (as compared with the one 
used for Moo — 1-55 and 2.0) were concentrated in regions where rapid changes of 
flow were expected like regions near shock waves, and the rear of the external tank 
(the sting was removed from the ET). Figures 4.3d and 4.3e show the surface Cp 
comparisons for the orbiter at the <j> = 70° line, and as indicated in the figures, the 
computations did not use the correct elevon deflection. That is the possible cause 
for the disagreement in the Cp comparisons at the trailing edge of the orbiter wing. 
Other than the orbiter wing trailing edge region, the numerical solutions are generally 
in better agreement with the wind tunnel data than with the flight data. 

At Moo = 1-05, a refined grid of 771,033 points was used which includes the 
ET/ORB forward and aft attach grids in addition to the ET, SRB, and ORB grids. 
The Reynolds number, Re = 4.0 x 10 ^ /ft, was taken from the wind tunnel test, and 
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(a) Top view 



(b) Side view 

Figure 4.2: Comparison of computational and wind tunnel surface pressure coeffi- 

cient at Moo = 1-55, a = —6°, and Re = 3.2 x 10 ^ / ft (3% model) 
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(a) Moo — 2.0 and a = —4°; (computation: 
no elevon deflection; flight: 8°/ — 5° in- 

board/outboard elevon deflection, a = —2°) 



(b) Moo = 1-55 and a = -6°, Re = 3.2 x 10 6 / ft; 

(computation: no elevon deflection; wind tun- 
nel: 10°/ — 7° elevon deflection; flight: 8°/ — 

5° inboard/outboard elevon deflection, a = 

- 2 °). 

Figure 4.3: Comparison of Cp from computation (-), wind tunnel (o), and flight ( V 

right side, A left side) along the <f> = 70° line of the orbiter fuselage at 
Moo = 2.0, 1.55, 1.05, 0.9 and 0.6 (x in inches) 
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(c) Moo — 105 and a — — 3°, (d) Moo = 0.9 and a = —3°, Re = 4.0 x 

Re = 4.0 x 10®//i; (all: 10°/9° in- 10®//*; (computation: no elevon 

board/outboard elevon deflection). deflection; wind tunnel and flight: 

10° /9° inboard/outboard elevon de- 
flection). 



10^//*; (computation: no elevon 

deflection; wind tunnel and flight: 
10° /9° inboard/outboard elevon de- 
flection). 


Figure 4.3 (Continued) 
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the angle of attack, a = -3°, and the elevon deflection, 10° /9° (inboard/outboard), 
were those of the wind tunnel and flight data. As observed from the surface pressure 
coefficients for the orbiter along the fuselage <f> = 70° line in Figure 4.3c, the Cp at 
the orbiter trailing edge is in better agreement than those of other Mach numbers. 
Overall, the computation is in better agreement with the flight data than with the 
wind tunnel data for most part of the fuselage — most likely due to the wall inter- 
ference in the wind tunnel test. Additional surface pressure comparisons between 
the numerical solutions and the wind tunnel data are presented in Figures 4.4 and 
4.5. The shaded surface pressure coefficient comparisons in Figure 4.4 show a simi- 
lar variation in the pressure contours for both the computational and experimental 
results. The quantitative comparisons of Cp at various constant angle lines for the 
ET, SRB and ORB are shown in Figure 4.5. For the most part, the computation is 
in good agreement with the wind tunnel data. The discrepancy near the back of the 
orbiter and the external tank, as observed in Figures 4.5a and 4.5b, can in part be 
attributed to the inadequate modeling of the ET/ORB attach hardware which only 
accounts for 50% of the blockage incurred from the real attach hardware and the 
fuel feed lines. The other cause for the discrepancy is that the wind tunnel model 
only had stings attached to the SRB nozzles while the numerical model had a sting 
extending from the back of the orbiter in addition to the SRB stings. For the Cp of 
the SRB in Figure 4.5c, the absence of the ET/SRB attach ring in the SRB grid is 
thought to be the major reason for the discrepancy between the numerical solutions 
and the wind tunnel data. 
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(a) Top view 



(b) Side view 

Figure 4.4: Comparison of computational and wind tunnel surface pressure coeffi- 

cient at Moo = 1-05, a = -3°, and Re = 4.0 x 10 6 //< (3% model) 
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(1) <f> = 0° (bottom) (2) <j> = 90° (side) 



(3) <t> = 180° (top) 


Figure 4.6: Comparison of C p from computation (-) and wind tunnel (o) for var- 

ious lines along the external tank for Moo = 1-05, a = —3°, and 
Re = 4.0 x 10 ®//< (x in inches) 
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(3) <f> — 180° (top) (4) $ — 270° (inner side) 

Figure 4.7: Comparison of C p from computation (-) and wind tunnel (o) for various 

lines along the solid rocket booster for Moo = 1-05, a = -3°, and 
Re = 4.0 x 10 ^ l ft (x in inches) 
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4.2 Vorticity- Velocity Solutions 

Numerical solutions were obtained for several geometries including a sphere, an 
ellipsoid, the ET, and the integrated space shuttle vehicle and compared with ex- 
perimental results, exact solutions, and numerical solutions of the thin-layer Navier- 
Stokes and Euler equations. Like the grids used for the thin-layer Navier-Stokes 
calculation, the grids used for the vorticity-velocity calculation were all generated 
using a hyperbolic grid generator (Steger and Rizk 1985). The sphere test cases are 
presented first to illustrate various aspects of the numerical scheme, such as the ac- 
curacy of the vorticity-velocity formulation, grid refinement effects, and effects of the 
outer boundary location on the solution. Then the ellipsoid testcases are presented 
for further validation of the numerical scheme. A comparison with experimental re- 
sults is included for the ellipsoid at 10° angle of attack. The solution of the ET alone 
is then presented followed by the solution of the integrated space shuttle vehicle. 

4.2.1 Sphere 

Since the analytical solution for a potential flow around a sphere is available 
for comparison, and the sphere grid is relatively easy to generate, the flow past a 
sphere was chosen as the first test case to verify the vorticity-velocity scheme. The 
flow was assumed inviscid and irrotational and the calculations were carried out for 
both the single as well as the chimera overset grids to verify the vorticity-velocity 
formulation in the chimera approach. Grids of different sizes (or different number 
of points) were used to investigate the effect of grid refinement. Since the far-field 
boundary condition used for the sphere is of the Dirichlet type (i.e., values of the flow 
variables are specified at the far-field boundary — in this case, free stream values 
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were used), it is important to set the far-field boundary far enough away from the 
surface to obtain an accurate solution. Thus, the effect of the far-field boundary 
location on the solution was evaluated. Comparisons of the Poisson solution with the 
analytical (Kaplan, 1940) and Euler solutions were made to evaluate the accuracy 
of the scheme. Flow conditions were set at Mach number 0.57 and zero angle of 
attack for a unit sphere. The symmetry plane boundary condition was used to save 
computer time. The tangency/vorticity consistency boundary condition was imposed 
on the body surface and the axis averaging boundary condition was used for the two 
axes extending from both poles of the sphere to the far-field boundary. For the Euler 
solution, the tangency boundary condition (explained previously in the chapter on 
“Governing Equations and Numerical Methods”) was used on the surface and the 
rest of the boundary conditions were the same as those used for the vorticity- velocity 
algorithm. 

The analytical solutions of a compressible, irrotational flow past a circular cylin- 
der and a sphere were first calculated by Janzen (1913) and Rayleigh (1916). Their 
method added a correction term, which involved only the square of the Mach num- 
ber, to the incompressible flow solution. Kaplan (1938) and Imai (1938) extended the 
calculations by including the terms involving the fourth power of the Mach number 
for a circular cylinder. 

<t> = <j> 0 + <t> 1 M 2 + <£ 2 M 4 + --- (4.1) 

2 

where <f> is the velocity potential and <j> q the solution of the Laplace equation, V 4> g = 

0, for an incompressible flow. By inserting Eq. 4.1 into the continuity equation, 

d^4> d^4> d^<j> 1 / dv ^ d<f> dv ^ d<j> dv ^ d<f>\ 

dx 2 + dp + dz 2 2c 2 \ dz dx dy dy + dz dz J 


(4.2) 
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(where v is the fluid velocity and can be expressed in terms of velocity potential, <j>) 
and equating the coefficients of the same powers of Mach number on both sides of 
the equation, <f> 0 , * f, and <f> 2 can be determined. Kaplan (1940) used this method 
to solve for the irrotational, compressible flow past a sphere. His solution for the 
velocities on the surface of the sphere is 

-2^ = | sin0 + (989 sin# - 1215 sin30)M 2 

+ (0.10572 sin 9 — 0.16008 sin 30 + 0.06434 sin50)M 4 ( 4 - 3 ) 

+ ( 7 - 1)(0. 01168 sin 9 — 0.02475 sin 30 + 0.02582 sin50)M 4 

where 7 is the ratio of the specific heats. For a flow at Mach number of 0.57, and 
7 = 1.408, the surface velocities on a sphere are: 

— = 1.55731 sin 9 — 0.07404 sin 30 + 0.00791 sin 50 (4.4) 

Voo 

4. 2. 1.1 Single grid calculation As illustrated in Figure 4.8, the sphere 
grid is clustered near the surface and quickly stretched out to the far-field boundary. 
The reason for using a viscous grid (grids with very fine grid spacings near the wall 
boundary) for the inviscid flow solver used here is to test whether the current scheme 
would work on such a grid since one of the major objectives for developing the 
vorticity- velocity formulation is to see if its solution can be used to provide better 
initial guess for a viscous flow solver such as one using the Navier-Stokes equations. 
If different grids are used for the vorticity-velocity algorithm and the viscous flow 
solver, say for the thin-layer Navier-Stokes equations, the error introduced from the 
interpolation could be large, especially in the boundary layer where most of the grid 
points reside. For the chimera overset grid with holes embedded, the interpolation 
will be much more complicated since the hole points contain no meaningful data and 
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Figure 4.8: Geometry of the computational domain for a sphere 
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the interpolation has to exclude such data from being used. Thus it is desirable to 
use the same grid for both the inviscid and viscous flow solvers. 

As indicated in Table 4.1, the flow around a sphere was computed for various 
far-field boundary locations, from 4 to 16 diameters from the center of the sphere and 
stretching ratios^- around 14.5%. The comparison of surface velocities for these cases 
are shown in Figure 4.9. There appears no appreciable difference in the solutions 
for different far-field boundary locations. Thus, for other calculations for the sphere 
geometry, the far-field boundary location was set at 4 diameters from the center of 
the sphere. Also shown in the same figure is the analytical solution of Kaplan (1940). 
The Poisson solution predicted lower velocities than that of Kaplan in regions near 
<f> = 90° (0.38% lower at 90°); however, in general, the two solutions agree well with 
each other. 

Table 4.1: Far-field boundary locations used for flows past a sphere 


Grid # 

Grid size (£ x 77 x () 

Far-field boundary 

Stretching ratio 

1 

39 x 39 x 47 

4 dia. 

14.5% 

2 

39 x 39 x 51 

6 dia. 

14.6% 

3 

39 x 39 x 53 

8 dia. 

14.4% 

4 

39 x 39 x 58 

16 dia. 

14.5% 


9 

Table 4.2 lists the size of several grids of different initial spacings ( ), ranging 
from 0.0005 to 0.025. The surface velocity from solutions of these five grids are shown 

*The “stretching ratio” listed in Table 4.1 in percentage actually represents the 
percentage increase in grid spacing from the body outward (the £ direction). That 
is, for a table entry of 14.5%, the true ratio would be 1.145. 

^A£j is the physical distance in the £ direction between the surface and the first 
point out in the flowfield normalized by the diameter of the sphere. 
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along with the analytical solution (Kaplan, 1940) in Figure 4.10. The solution on 
the coarsest grid (21 x 21 x 24) is in total disagreement with the analytical solution. 
With the second coarsest grid (27 x 27 x 30), the surface velocities are much closer 
to the analytical solution, but still somewhat larger than that given by analytical 
solution. On an even finer grid (33 x 33 x 35) with an even smaller initial spacing 
(A(i = 0.005), the solution generally agrees with the analytical solution except in 
regions near 90° where the Poisson solution predicts somewhat larger velocities than 
those of the analytical solution. The solutions from the two finest grids (39 x 39 x 47 
and 47 x 47 x 52) are almost identical suggesting that the solution on the finest grid 
can be considered as the “true” solution of the “finite-difference Poisson equations”. 
Although the solutions from the two finest grids generally under-predict slightly the 
surface velocity compared to the analytical solution, they are in reasonably good 
agreement with the analytical solution. 


Table 4.2: Grid sizes used for flows past a sphere 


Grid # 

Grid size 

u 

X 

V x C) A<1 

Stretching ratio 

1 

21 

X 

21 

X 

24 

0.025 

14.1% 

2 

27 

X 

27 

X 

30 

0.01 

14.6% 

3 

33 

X 

33 

X 

35 

0.005 

14.6% 

4 

39 

X 

39 

X 

47 

0.001 

14.5% 

5 

47 

X 

47 

X 

52 

0.0005 

14.5% 


An Euler solution was computed for a grid of 39(4) x 39(t/) x 47(4) points with 
the initial spacing, = 0.001, in the radial direction and a 14.5% stretching ratio. 
Figure 4.11 shows surface velocities from solutions of the Poisson and Euler equations 
as well as those from the analytical solution for compressible (Kaplan, 1940) and 
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Figure 4.9: The effect of the far-field boundary location on the Poisson solution for a sphere at Moo 
and a = 0° 
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Figure 4.10: The effect of grid refinement on the Poisson solution for a sphere at Moo = 0.57, and a = 0° 
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Figure 4.11: Comparison of the Poisson and Euler solutions with the analytical solution for a sphere at Moo = 
0.57, and a = 0° 
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incompressible flows. Overall, the Poisson and Euler solutions agree well with the 
analytical solution for the compressible flow (Kaplan, 1940). The difference between 
the Poisson and the Euler solutions is that the Poisson solution under-predicts while 
the Euler solution over-predicts the analytical solution for the compressible flow near 
the 90° region. Also observed for the Euler solution is the slight asymmetry in the 
solution as seen from the minor downstream “shift” from Kaplan’s analytical solution. 
The major discrepancy between the solutions of the compressible and incompressible 
flows is in regions near 90° where as much as 8.3% difference is observed. 

Figure 4.12 shows the typical convergence history of the Poisson solver imple- 
mented in the present study. The L2-norm is found to drop very quickly in the early 
stages of the computation; then the rate of convergence slows down. However, the 
L2-norm usually drops two orders of magnitude in 100 iterations. The CPU time 
required for the Poisson solver is about 4.67 fis per point per iteration on a Cray 
YMP and requires about 150 to 500 iterations, depending on the size of the grid, 
to obtain a converged solution. The Euler solver requires about 31.58 n s per point 
per iteration and requires about 1000 iterations to converge for the same grid. The 
convergence criterion is based on the magnitude of the L2-norm. If it is three orders 
of magnitude smaller than the peak value during the solution process, the solution is 
considered converged. 

4. 2. 1.2 Chimera grid calculation For the chimera scheme, two overset 
grids were computed and compared with the single grid solution. Both the chimera 
and single grid solutions were expected to be the same or very close to each other. 
The flow was maintained at a Mach number 0.57, and zero angle of attack. All 
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Figure 4.12: Convergence history of the Poisson solution for a sphere at Moo — 0.57, and a = 0° 
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Figure 4.13: Surface velocity comparisons of the single and chimera grid solutions for a sphere at Moo = 0.57 
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grids (chimera or single) used were uniformly spaced in the £ (streamwise) and 77 
(circumferential) directions and stretched to the far-field boundary (4 diameters from 
the center of the sphere) from an initial spacing of 0.001 diameters in the ( (radial) 
direction. The single grid had 39, 39, and 47 points in the £, 77 and < directions, 
respectively. The two overset grids were taken from subsets of the single grid to 
avoid the possible error introduced from the interpolation at the chimera interface 
boundaries. The first one, designated as chimera grid I in Figure 4.13, consisted of 
the first thirty (£ = 1 to 30) and the last twenty-one (£ = 27 to 47) constant ( planes 
of the single grid while the second one, designated as chimera grid II, consisted of the 
first twenty-two (£ = 1 to 22) and the last twenty-two (£ = 18 to 39) constant £ planes. 
The boundary conditions for the chimera grid were the same as those for the single 
grid testcase, i.e., free stream values at the far-field boundary, tangency/vorticity 
consistency boundary conditions on the surface of the sphere, and the axis averaging 
condition for the axes. Comparison of surface velocities from the single and chimera 
grids are shown in Figure 4.13. For the chimera grid I (four points overlapped in 
the ( direction), the surface velocity at 90° is 1.57% lower than that of the single 
grid. This discrepancy is due to the error resulting from the one-sided differencing 
used in computing the metrics at the chimera interface boundary compounded with 
the grid stretching in the < direction. While for the chimera grid II, the one-sided 
differencing used in calculating the metrics at the chimera interface boundaries does 
not affect the accuracy since the chimera interface boundaries are uniformly spaced 
£ planes. This is evidenced by the very good agreement between the solutions of the 
chimera grid II and the single grid. The minor disagreement (0.3% difference) at 90 
is likely due to the different differencing scheme used in the chimera scheme to take 
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into account those points at the chimera interface boundary. 

For the cases computed, the chimera scheme required about 4.5 times the CPU 
time (20.94 jts vs. 4.67 ps per point per iteration on a Cray YMP) needed for a 
corresponding single grid solution due to the extra logic required to identify the hole 
and the chimera interface boundary points from the field points and the extra disk 1/ 0 
to bring in and put back the data from/to the disk. The time for disk I/O represented 
30% of the total CPU time used and could have been reduced to a bare minimum if 
the SSD (solid state device) was used. Thus the actual CPU time required per point 
per iteration, not counting the disk I/O, for the chimera scheme is about 3.1 times 
that of a single grid. However, due to the explicit type boundary conditions used at 
the chimera interface boundary, an additional 20% more iterations were required to 
converge the solution than for the corresponding single grid solution. 

4.2.2 Ellipsoid 

The geometry of the computational domain for the ellipsoid computed is shown in 
Figure 4.14; the grid is clustered near the wall and the first grid spacing away from the 
surface is one ten-thousandth of the length of the ellipsoid. Unless stated otherwise, 
all the cases computed for the ellipsoid in this research used a 39(C) x 21 (t/) x 51(C) 
grid. The £ direction is aligned with the incoming flow, q, the circumferential di- 
rection, and C) the direction away from the body surface. The aspect ratio of the 
ellipsoid is 6 to 1 (length vs. diameter of the ellipsoid), and the far-field boundary, a 
near spherical surface, is at a distance about 5.4 times the length of the ellipsoid from 
the center. A symmetry plane condition is used to save CPU time. The wall bound- 
ary condition is the usual no-slip boundary condition for the thin-layer Navier-Stokes 
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solutions, the tangency (no penetration) boundary condition for the Euler solutions, 
and the tangency /vorticity consistency boundary condition for the Poisson solutions. 
The axis boundary condition is used for the two axes extending from both ends of 
the ellipsoid to the far-held boundary. For details on these boundary conditions, see 
Chapter 3 on the “Governing Equations and Numerical Methods”. 

In the 0° and —3° angles of attack testcases, the Mach number was taken as 
0.6 and the Reynolds number, 10 million. For the thin-layer Navier-Stokes solution, 
the flow was assumed turbulent and the Baldwin- Lomax turbulence model was used. 
Figure 4.15 compares the surface pressure coefficients from the thin-layer Navier- 
Stokes, Euler and Poisson solutions for an ellipsoid at the plane of symmetry for 
a flow at 0° angle of attack. All three solutions show good agreement with each 
other in most of the flow region except far downstream where the flow separates due 
to an adverse pressure gradient. The pressure coefficients on the top and bottom 
surfaces of the ellipsoid are found to be the same for all three solutions as expected 
for axisymmetric flows. A similar comparison is also presented in Figure 4.16 for a 
plane perpendicular (<j> — 90°) to the plane of symmetry. The solution at this plane 
is essentially the same as the solution on the plane of symmetry due to 0° angle of 
attack, and thus the same level of agreement is observed. For the a = — 3° case, the 
thin-layer Navier-Stokes solutions were computed and compared in Figure 4.17 for 
the symmetry plane, and in Figure 4.18 for the 4> = 90° plane. The Poisson, Euler, 
and thin-layer Navier-Stokes solutions agree with each other in most regions of the 
flowfield although the Poisson solution gives somewhat lower pressure than the other 
two solutions. The major difference is at the downstream end of the ellipsoid due to 
the separation resulting from the adverse pressure gradient. Unlike the solution for 


Figure 4.14: Geometry of the computational domain for an ellipsoid 
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the 0° angle of attack case, the Poisson solution in this case is not symmetric due to 
the non-zero angle of attack. 

The third case for the ellipsoid geometry is intended to validate the vorticity- 
velocity algorithm with the experiment done by Kreplin, Vollmers and Meier (1982), 
and Meier and Cebeci (1985). The Reynolds number based on the free stream condi- 
tions and the major axis is 1.6 x 10 6 , the Mach number, 0.4, and the angle of attack, 
10°. Since the experimental study indicated that the flow was nominally laminar, the 
simulation was carried out for laminar flow. However, unlike the geometry used in 
the numerical simulation, the real experimental setup had a sting attached to the end 
of the ellipsoid. Vatsa, Thomas and Wedan (1987) computed a Navier-Stokes solu- 
tion for this case using upwind-biased and central differencing for the convective and 
pressure terms respectively, with and without the sting for laminar and transitional 
flows. Their study showed that the sting mostly affected the flow near the juncture of 
the sting, and the flow features for the rest of the flowfield were nearly the same with 
and without the sting. Thus, the current calculation was carried out without the 
sting. From the comparison of the surface pressure coefficients in Figure 4.19, it can 
be observed that the Poisson solution generally predicts a lower pressure than was 
found in Vatsa, Thomas and Wedan (1987) for both the windward and the leeward 
sides. Nonetheless, the Poisson solution agrees well with their results for most of 
the leeward side and the first half of the windward side except near the downstream 
end of the ellipsoid where the flow separated from the surface as illustrated by the 
particle traces (Figure 4.20) for the thin-layer Navier-Stokes solution. The thin-layer 
Navier-Stokes solution, like the Poisson solution, predicts lower pressure than the 
experimental results and is in good agreement with the experiments on the first half 
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Figure 4.16: Comparison of Cp from the Poisson, Euler and Navier-Stokes solutions at the <j> = 90° plane for an 
ellipsoid at Moo = 0.6, Re = 10^, and a — 0° 












Cp for 6:1 EllLpsold at k = 1 1 

Rlpha-~3, Mach-0.6, Re~10000000, Grld _ 39x21x51 


LEGEND 

NS 

Euler 


PoLsson 


Figure 4.18: Comparison of Cp from the Poisson, Euler and Navier-Stokes solutions at the <f> = 90° plane for an 
ellipsoid at Moo = 0.6, Re = 10^, and a = — 3° 
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of the windward side and disagrees at the downstream end where the computational 
model differs from the experimental setup. However, the flow characteristics for the 
thin-layer Navier-Stokes solution at the downstream end is similar to that of Vatsa, 
Thomas and Wedan (1987) where the leeward pressure was higher than the wind- 
ward pressure. Overall, the thin-layer Navier-Stokes solution is quantitatively similar 
(though not perfect) to the experimental results and the level of agreement can be 
improved by increasing number of grid points and reducing the grid spacing in the 
boundary layer to resolve the gradients. In Vatsa, Thomas and Wedan (1987), the 
grids used were 4.2 times (73 x 49 x 49) for the case with the sting and 11 times 
(97 x 97 x 49) for the sting-less case as many points as used in the current study. 
Since the primary purpose of the numerical calculation for this case was to demon- 
strate the capability of the vorticity-velocity algorithm, the thin-layer Navier-Stokes 
solution was not pursued further to obtain better agreement with the experimental 
results. 

One of the objectives of the current research was to find ways to improve the effi- 
ciency of the numerical methods for flow simulations. The proposed method is to first 
obtain the flow solution from a set of simplified equations (here, the vorticity-velocity 
formulation is used with vorticity being set to zero, i.e., the Poisson equations) and 
then either directly feed that solution into the Navier-Stokes solver or “combine” the 
Poisson solution with the Navier-Stokes equations to obtain an improved solution 
and then solve the Navier-Stokes equations for the entire flowfield. To combine the 
Poisson solution with the Navier-Stokes equations means that the Navier-Stokes (or 
the thin-layer Navier-Stokes) equations are solved for a user specified region which 
roughly encompasses the boundary layer, using the Poisson solution as the edge 
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Figure 4.19: Comparison of Cp from the experimental results, the Poisson and Navier-Stokes solutions at the 
plane of symmetry for an ellipsoid at M 0 0 — 0.4, Re — 1.6 x 10^, and a — 10° 
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Figure 4.20: 


Particle traces from the thin-layer Navier-Stokes solutions for an ellip- 
soid at Mo o = 0.4, Re = 1.6 x 10 6 , and a = 10° 
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(c) Leeward, thin-layer Navier-Stokes solutions. 
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(d) Leeward, combined Poisson and thin-layer Navier-Stokes 
solutions. 


Figure 4.21 (Continued) 
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Figure 4.22 (Continued) 
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conditions. This step need not be. carried out until convergence since it is only an 
intermediate step in obtaining an even better estimate than provided by the Pois- 
son solutions alone. However, from the cases tested, it was found that solving the 
Navier-Stokes equations for half (or more) of the grid points, the rate of convergence 
(based on the residual) was an order of magnitude faster than solving the Navier- 
Stokes equations for the entire flowfield. This clearly shows that the deterioration 
of convergence is faster than the linear increase of number of grid points used and 
the larger the grid, the more CPU time will likely be saved by improving the initial 

guess. 

In Figure 4.21, the convergence history of the pressure coefficient is shown for the 
thin- layer Navier-Stokes solutions using a uniform flow field as the initial guess and 
also using the combined Poisson and thin-layer Navier-Stokes solutions as the initial 
guess before solving the entire flowfield using the thin-layer Navier-Stokes equations 
for the Moo = 0.4, Re = 1.6 x 10 6 and a = 10° testcase. Use of the better initial guess 
appears to significantly reduce the oscillations in the solution process compared to the 
thin-layer Navier-Stokes results obtained from using a uniform flowfield as the initial 
guess. Note that the uniform flowfield referred to here was actually implemented with 
a “slow-start” mechanism which gradually reduced the velocities on the surface from 
the free stream values to zero to avoid the possible instability resulting from suddenly 
setting the wall velocities to zero. From here on, the uniform flowfield referred to 
for the initial guess implies that the “slow start” is used. For both the windward 
and the leeward sides, the better initial guess required 1073 iterations to converge 
while the one using the uniform flowfield as the initial guess needed 1790 iterations 
to achieve the same level of convergence. Adding the CPU time required to obtain 
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the better initial guess, which is equivalent to 145 iterations of the thin-layer Navier- 
Stokes solutions, the total savings of CPU time is about 32% (or 572 iterations). For 
the other case, Moo — 0.6, Re = 10^ and a = —3°, the Poisson solution (without 
using the combined scheme) was used as the initial guess to obtain the thin-layer 
Navier-Stokes solution and is compared with the thin-layer Navier-Stokes solution 
computed initially from a uniform flowfield in Figure 4.22. Overall, when the Poisson 
solution was used as the initial guess, smaller oscillations in the solution process were 
observed than when the a uniform flowfield was used as the initial guess. With the 
better initial guess, a converged solution was obtained in 1390 iterations as shown in 
Figure 4.22 while 2290 iterations was required when the uniform flowfield was used 
as the initial guess. Adding the CPU time required to obtain the Poisson solution, 
which is equivalent to 82 iterations of the thin-layer Navier-Stokes solutions, the 
savings of CPU time is 36%. All the cases presented here were computed on a Cray 
YMP computer, and the CPU time required for each point per iteration is 43.8 fis 
using the thin-layer Navier-Stokes equations and 9.98 \i s for the Poisson equations. 

4.2.3 External Tank (ET) 

Since the external tank is one of the major components of the integrated space 
shuttle vehicle, and its grid was used as the major grid to cover the entire compu- 
tational domain in the shuttle flow simulation, it was used as yet another testcase 
before delving into the shuttle flow simulation using the vorticity- velocity algorithm. 
Two different angles of attack, 0° and -3°, were calculated for a flow of M 0 o = 0.2, 
and Re = 4.0 x 10°//t. Both the thin-layer Navier-Stokes and Poisson solutions were 
carried out and compared. As shown in Figure 4.23, the computational domain of the 
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external tank grid starts from three times the body length upstream to 2.7 times the 
body length downstream of the ET and the far-field boundary is at least three times 
the body length from the ET center of gravity. Also noticed in Figure 4.23 is a sting 
approximately the diameter of the ET extending from the back end of the ET to the 
downstream boundary. The sting was used to circumvent the inability of the Poisson 
equation to model the separated flow. The grid used 63, 39, and 51 points in the £ 
(streamwise), 77 (circumferential), and C (outward-going) directions, respectively, and 
like the rest of the grids used in this study, most of the points were clustered near the 
body surface. The initial spacing away from the body surface was 5.936 x 10“^ times 
the ET length (0.13/t). The grid points clustered near the rear end of the external 
tank in the £ direction were intended to resolve the flow around the aft attach hard- 
ware and were not actually needed here. However, they were used for convenience. 
The wall boundary conditions used were no-slip and tangency/vorticity consistency 
conditions for the thin-layer Navier-Stokes and the Poisson solutions, respectively. 
The axis boundary condition was used for the axis extending upstream from the ET 
nose to the far-field boundary. Free stream values of flow variables were used at the 
far-field boundary, and a zero-gradient outflow condition was used at the downstream 
boundary. A symmetry plane boundary condition was employed in the 77 direction 
to save CPU time. 

For the a = 0° case, the surface pressure coefficient comparisons are presented 
in Figure 4.24. The thin-layer Navier-Stokes solutions were computed with a uniform 
flowfield as the initial guess. Since this flow is axisymmetric due to zero angle of 
attack, the pressure at the top and bottom of the symmetry plane are the same, as 
shown in the figure. The thin-layer Navier-Stokes solutions show a sharper expansion 
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Figure 4.23: Geometry of the computational domain for the external tank (x in 10^ 

inches) 
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Figure 4.24: Comparison of surface pressure coefficients for the external tank, Moo = 0-2, Re = 4.0 x 10 °//t 
and a = 0° (x in 10^ inches) 




Cp for ET ( single grid) at plane of symmetry 


Alpha- 3, Mach-0 . 2 , Re— 10000000, Grld-63x39x51 




09 




0 1 2 3 4 5 6 7 

X 


(b) Windward, combined Poisson and thin-layer Navier- 
Stokes solutions, a = —3° 


Figure 4.26: Convergence history for the external tank at M 0 © = 0.2, and 

Re = 4.0 x 10 ^ /ft at iteration 190 ( ), 590 ( ), 1190 ( ), 

1790 (* • •), 2390 ( — ), (x in 10^ inches) 
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(c) Leeward, combined Poisson and thin-layer Navier-Stokes 
solutions, a = — 3° 

Figure 4.26 (Continued) 
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than the Poisson solutions at the front portion of the external tank followed by a 
wiggle in the solution which is likely due to the lack of grid points to resolve the 
flow. Then both solutions return back to free stream pressure for the rest of the 
external tank. At the juncture of the external tank and the sting, the pressure from 
the Poisson solution shows a sharper variation than that of the thin-layer Navier- 
Stokes solutions due to lack of viscous terms to smooth out the solutions. Then, 
both solutions gradually return to the free stream pressure for the rest of the sting. 
In Figure 4.25, the same trend — sharper expansion for the thin-layer Navier-Stokes 
solutions in the front portion of the ET and the sharper pressure variation for the 
Poisson solutions at the juncture of the ET and the sting — can be found for the 
a = —3° case. However, the top and bottom sides of the symmetry plane have 
different pressures due to non-zero angle of attack. Generally speaking, for both 
cases presented here, the Poisson solutions agree with the thin-layer Navier-Stokes 
solutions for most regions of the external tank and can be considered as a good initial 
guess for the thin-layer Navier-Stokes solutions if such a solution is desired. 


The convergence history for the external tank is shown in Figure 4.26. In Fig- 
ure 4.26a, the convergence history of the thin-layer Navier-Stokes solution computed 
from a uniform flow field with a slow start shows large oscillations during the solu- 
tion process and is not converged yet after 2390 iterations while for the case using 
the combined Poisson and the thin-layer Navier-Stokes solutions as the initial guess, 
shown in Figures 4.26b and 4.26c, the level of oscillation is significantly smaller and 
the solution converged in 1790 iterations. 
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4.2.4 Integrated Space Shuttle Vehicle (ISSV) 

The primary objective of developing the vorticity-velocity algorithm was to be 
able to economically simulate the flow around the integrated space shuttle vehicle. 
Since the algorithm essentially involves solving the Poisson equations, it is faster than 
solving the Navier-Stokes or the thin-layer Navier-Stokes equations, and thus can be 
used to evaluate the overset grid used in the chimera scheme with only modest use 
of computing resources. The resulting solution can also be used as an initial guess 
for more accurate algorithms which solve the Navier-Stokes or the thin-layer Navier- 
Stokes equations. 

Due to limited computing resources, only the accuracy of the vorticity-velocity 
algorithm is evaluated in the chimera approach. The convergence acceleration using 
the Poisson solution as the initial guess was not carried out; however, it is believed 
that the rate of convergence will be improved using the Poisson solution as the initial 
guess as was found for the ellipsoid and external tank geometries. Like the grids 
used earlier for the thin-layer Navier-Stokes flow simulation of the integrated space 
shuttle vehicle, the ET grid was made the major grid which extended all the way 
out to the far-field boundary while the SRB and ORB grids occupied a much smaller 
portion of the computational domain. The computational domain started roughly 
3 times the ET length upstream to 2.7 times the ET length downstream of the ET 
with a grid of 383,001 points spread through the three component grids. All three 
grids were generated in the same manner as described in the section on the thin-layer 
Navier-Stokes solutions in this chapter. Since the purpose of computing the integrated 
space shuttle vehicle here was to verify the vorticity-velocity algorithm in the chimera 
approach, the grids used were not intended to reflect the true space shuttle geometry. 
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For example, the ET/ORB attach hardware was not modeled in the geometry and 
all three component grids, ET, SRB and ORB, had stings (including the orbiter wing 
sting) extending from the back of each component to the downstream boundary. 

As with the thin-layer Navier-Stokes equations, uniform free stream flow was 
assumed initially for the whole computational domain in the Poisson solution pro- 
cedure. The solution process started first by solving flows on the ET grid with the 
proper boundary conditions. The flow variables on the hole boundary were still at 
free stream values at this stage. After the ET solution was obtained, the variables 
on the outer and hole boundaries of the SRB and ORB grids were then updated 
by the ET solution. The SRB grid was then solved with the newly updated outer 
and hole boundary values as well as other necessary boundary conditions. Then this 
SRB solution was used to update the flow variables on the hole boundaries of the 
ET and ORB grids and the outer boundary of the ORB grid. Likewise, the solution 
on the ORB grid was obtained and used to update the other two grids. This process 
repeated until convergence was reached. Although it is not mentioned above, only 
one grid was solved at a time and the solution was put back to external disk or SSD 
before another grid was brought into the computer core memory and solved. Unlike 
the thin-layer Navier-Stokes solution procedure, the slow start was not required for 
the Poisson solution procedure since the wall velocities did not suffer from an abrupt 
change of values as in viscous flows. 

Both the thin-layer Navier-Stokes and Poisson solutions were carried out and 
compared. Free stream values were assumed on the ET far-field boundary. The 
downstream boundaries of all three grids used a zero-gradient boundary condition 
which forced the flow variables at the last two downstream stations to be the same. 
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Figure 4.27: Comparison of surface pressure coefficients from the Poisson and 

thin-layer Navier-Stokes solutions for the ET in ISSV configuration 
(x in 10^ inches) 
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(a) 0 = 0° (bottom, leeward) 
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(b) <fr — 180° (top, windward) 

Figure 4.28: Comparison of surface pressure coefficients from the Poisson and 

thin-layer Navier-Stokes solutions for the SRB in ISSV configuration 
(x in 103 inches) 
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Figure 4.28 (Continued) 
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Figure 4.29: Comparison of surface pressure coefficients from the Poisson and 

thin-layer Navier-Stokes solutions for the ORB in ISSV configuration 
(x in 10 3 inches) 
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Cp for ORB (sLde) 
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Figure 4.29 (Continued) 




120 


The values on the axes extending from the ET, SRB and ORB noses to the upstream 
far-held boundary were calculated by taking the average from points surrounding 
the axes. On the body surface, the tangency/vorticity consistency and the no-slip 
boundary conditions were employed for the Poisson and thin-layer Navier-Stokes 
solutions, respectively. The flow was computed for a Mach number of 0.6, and an 
angle of attack of —3°. Surface pressure coefficient comparisons for the ET, SRB 
and ORB at various constant angle lines are presented in Figures 4.27 to 4.29. For 
the ET, the thin-layer Navier-Stokes solution predicted a sharper expansion than the 
Poisson solution near the nose region; however, the two solutions are in a reasonable 
agreement with each other. At <f> = 180°, x = 1.25, and <f> = 90°, x = 1.0, the 
thin-layer Navier-Stokes solution also shows a sharper expansion than the Poisson 
solution, due to the sharper expansion predicted by the thin-layer Navier-Stokes 
solution at the orbiter nose and the SRB nose, respectively. At <f> = 180° and 
x R 2 2.0, the Poisson solution predicted higher pressures than the thin-layer Navier- 
Stokes solution due to the higher pressure predicted by the Poisson equations on the 
orbiter surface at the corresponding location. At the rear end of the ET (x % 2.4) 
for all three constant angle lines ( <j> = 0, 90, and 180°), the Poisson solution again 
shows higher pressure than the thin-layer Navier-Stokes solution. This is partly 
due to the higher pressure predicted by the Poisson equations at the OMS (orbiter 
maneuvering system) pod of the orbiter and near the nozzle of the SRB, and partly 
due to the Poisson solution trying to return to the stagnation pressure due to the 
shrinkage at the rear of the ET. For the SRB and ORB, a sharper expansion is 
observed for the thin-layer Navier-Stokes solution near the nose regions while higher 
pressure was predicted by the Poisson equations at the OMS pod of the orbiter and 
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near the nozzle of the SRB ( x sa 2.4). For the rest of the flow, the Poisson solution 
predicted higher pressure than the thin-layer Navier-Stokes solution. Generally, the 
farther away from regions where two components are near each other, the better is 
the agreement observed between the Poisson and thin-layer Navier-Stokes solutions, 
e.g., at <f> = 0° (bottom) of the ET, <f> = 180° (top) of the ORB and regions near 
noses of all three grids. This can be attributed to the flow inside the small clearance 
between the ET/SRB and ET/ORB being largely (or entirely) inside the boundary 
layer where the potential flow assumption is not valid. Overall, the surface pressure 
coefficients from the Poisson solution are quite different from those of the thin-layer 
Navier-Stokes solution. Thus, viscous effects need to be incorporated into the Poisson 
equations to capture meaningful flow features for the integrated space shuttle vehicle. 
This can be done through the viscous-inviscid interaction between the Poisson (or the 
vorticity- velocity) equations and the Navier-Stokes (or the boundary-layer) equations. 
This part, however, is not carried out in the current research and remains a future 
research topic. 

In Figure 4.30, the gray scale contours of the surface pressure coefficient show 
the three-dimensional effects of pressure variation for the two solutions. Like the 
quantitative comparisons of surface pressure coefficient at various constant angle 
lines, only regions near the noses of all three grids show similar patterns of pressure 
variations. As for the CPU time usage, the Poisson solution took 6850 seconds and 
960 iterations to converge while the thin-layer Navier-Stokes solution took more than 
16 hours of Cray-2 time and 2000 iterations to converge. 
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(a) Top view 

Figure 4.30: Comparison of C p from the Poisson and thin-layer Navier-Stokes solu- 

tions for ISSV at Moo = 0.6, Re = 4 x 10 6 //f, and a = -3° 



(b) Side view 


Figure 4.30 (Continued) 
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5. CONCLUSIONS 

Numerical simulations using the thin-layer Navier-Stokes equations and chimera 
approach for flows around the integrated space shuttle vehicle were carried out. In 
general, the computed surface pressure from the numerical solution was in good 
agreement with the available wind tunnel and flight data. However, there were dis- 
crepancies due to simplification or omission of geometric features, and improvements 
are needed for the numerical solution to be of greater use in engineering applications. 
Among the crucial improvements needed are the refinement in the modeling of the 
ET/ORB attach hardware and the addition of the ET/SRB attach ring. The ideal- 
ization or the lack of these parts in the numerical model is believed to have caused 
the major discrepancies between the numerical solutions and the wind tunnel and/or 
flight data. In the expansion region, either increasing the grid resolution or using 
a higher order (more accurate) scheme may improve the accuracy of the numerical 
solutions. Currently, a joint effort between NASA Johnson Space Center and NASA 
Ames Research Center to improve the computational model of the space shuttle ge- 
ometry is underway, e.g., the addition of the orbiter vertical tail and the space shuttle 
main engines. For the chimera scheme, using a conservative interpolation procedure 
instead of the current non-conservative trilinear interpolation is likely to help the 
shock capturing for the transonic or supersonic flows. 
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A Poisson solution procedure (a special case of the vorticity- velocity algorithm) 
using primitive variables was implemented to enhance the efficiency of the numerical 
flow simulation. Both single and chimera overset grids can be used with this Poisson 
solution procedure. The accuracy of the Poisson solution was validated^ with the 
analytical (Kaplan, 1940) and Euler solutions for a sphere geometry. The effects of 
the far-field boundary location and grid refinement were evaluated. It was found that 
a far-field boundary placed at 4 diameters from the center of the sphere was sufficient 
to obtain an accurate solution. The computed results were found to approach fixed 
values when the grid was refined. In the chimera approach, the solution for a sphere 
using chimera overset grids was found to be in good agreement with the corresponding 
single grid solution. However, the chimera overset grid solution required several times 
more CPU time than a corresponding single grid solution. This can be improved by 
reducing disk I/O time through the use of the SSD (solid state device) and optimizing 
the finite-differencing scheme in distinguishing the hole points from the field points. 

Additional validation of the Poisson solution was carried out for the ellipsoid 
and the external tank geometries by comparing it with the thin-layer Navier-Stokes 
and/or Euler solutions. In general, good agreement with the thin-layer Navier-Stokes 
and Euler solutions was observed for the Poisson solution. A Poisson solution for an 
ellipsoid at 10° angle of attack was computed and compared with the thin-layer 
Navier-Stokes solution as well as experimental results (Kreplin et al., 1982; Meier 
and Cebeci, 1985). Good agreement was observed for most regions except near the 
downstream end of the ellipsoid where flow separation occurred. By using the Poisson 
or the combined Poisson/thin-layer Navier-Stokes solutions as an initial guess for 

^Only surface pressure coefficients were compared. 
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the thin-layer Navier-Stokes equations, roughly 30% CPU time savings was found 
as compared with solving the thin-layer Navier-Stokes from a constant free stream 
flowfield. 

Finally, the Poisson equations were solved for the integrated space shuttle vehicle. 
The Poisson solution in general did not agree well with the thin-layer Navier-Stokes 
solution. This is believed to be due to the fact that the clearance between different 
component grids was very small so that the flows in these regions were mostly (if 
not entirely) inside the boundary layer, which the Poisson equations were incapable 
of modeling. Thus, some form of viscous-inviscid interaction (e.g., using the Poisson 
and thin-layer Navier-Stokes equations for viscous-inviscid interaction) is needed to 
capture the viscous effects in the boundary layer. A simpler way to do the viscous- 
inviscid interaction without computing the boundary-layer thickness is to solve the 
thin-layer Navier-Stokes equations for a user specified region and to solve the Poisson 
equations for the region that lies outside the layer in which the thin-layer Navier- 
Stokes equations are solved. The vorticity can be added to the Poisson equations to 
account for the rotational flow. 
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7. APPENDIX A. POISSON EQUATIONS IN GENERALIZED 

COORDINATES 


Poisson equations: 


V^u 4- $x + wgy — u>2 z = 0 

(7.1) 

V^v + ■dy 4- - uJ^X = ® 

(7.2) 

V^u> 4 i?z 4- ^2x — ^ly 0 

(7.3) 

s in vector form: 


d x E + dyF 4 dzG + H = 0 

(7.4) 


where 


,F = 




fix 4 ^3^ — w 2z 

H= dy + W lz - U» 3a . ( 7>6 ) 

t?z + W 2x , 

In generalized coordinates, the above equation can be written as 
(zE( +VxE n + CxE c +(yF^ + v y F v 4 (yF{ + tzG ^ 4- r,zG v 4- CzG( + H = 0 (7.7) 


BUI* 


PRECEDING PAGE BLANK NOT FILMED 


140 


Identities resulting from the coordinate transformation: 



Scale the Poisson equations by the Jacobian and combine them with the identities 
from coordinate transformation to form the Poisson equations in conservative form 


as shown below. 

+ JE V + <fB ( + E [(£) { + (2f), + ) c ] + 

( i F ( + V i F n + + F [$){ + {% + (^)J + < T - U > 

fa + Jo, ^G (+ G [(^ + (?), + (£) J + f = 0 

From the chain rule of differentiation, the above equation can be rewritten as, 


txE + tyF + ZzG\ fyxE + r) y F + tj z G\ ( ( x E + ( y F + ( z G\ H _ 

j )t\ j y/v j IS'-" 

( 7 . 12 ) 


where 




( 7 . 13 ) 


( 7 . 14 ) 
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8. APPENDIX B. IMPLICIT PROCEDURE FOR 
TANGENCY/VORTICITY CONSISTENCY BOUNDARY 
CONDITION WITH POISSON EQUATIONS 


The tangency/vorticity consistency boundary condition for the solid boundary 
is solved implicitly with the Poisson equations. The equation below illustrates that 
the finite-difference equations of the computational domain are a system of block 
tridiagonal matrices. Generally, a block tridiagonal solver is used to solve a block 
tridiagonal system of equations, and the number of operations required to invert this 
system of equations is approximately (5iVM^/3) where N is the number of equations 
and M the blocksize (Fletcher, 1988). However, due to the special structure of the 
matrix involved here, it is possible to apply the Thomas algorithm (used to invert 
scalar tridiagonal matrices) to invert the block tridiagonal system of equations. Since 
the number of operations needed to invert scalar tridiagonal systems is only about 
(5 NM) (Fletcher, 1988), the savings in number of operations is roughly (A/^/3 — 1). 
For the Poisson equation ( M = 3), so only one-third of the operations required to 
invert the block tridiagonal matrix is needed to invert the following block tridiagonal 
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matrices using the Thomas algorithm. 


Bi C x 
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( 8 . 1 ) 


^n— 1 ^n— 1 C n — 1 
An Bn 

In the above equation, A n , B n and C n are 3 by 3 diagonal matrices, except for C\ 
which is a full 3 by 3 matrix, and X n and D n are vectors of length 3, n = 1, 2, • • • , n. 
Since all matrices in the above equation are diagonal matrices (except Cj), it is easy 
to eliminate the upper diagonal matrices by Gauss elimination for n = n — 1, • • • , 2. 
Thus, by applying Gauss elimination, the above equation becomes 
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with the first rows of matrix elements assuming the form as illustrated below: 
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The solution of the above equation, and X 2 , can be obtained analytically. 
Then, a forward substitution can be used to get the solution for X n ,n = 3, • • • ,n. 



